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Foreword 


This  is  the  Ph.D.  thesis  of  Tommy  Jensen.  A  new  ocean  model  for  the 
MASIG  team  at  FSU  is  described.  It  is  applied  to  the  seasonal  circulation  of  the 
upper  Indian  Ocean.  The  model  is  a  4-layer  isopycnic  model  with  the  lowest  layer 
at  rest.  Such  a  model  produces  hundreds  of  millions  of  numbers.  It  is  not  easy 
to  understand  and  describe  the  entire  calculation.  Thus  we  decided  to  concentrate 
on  the  western  boundary  structure  in  the  Arabian  Sea  and  the  equatorial  seasonal 
currents. 

It  is  shown  that  a  reduced  gravity  model  isn’t  too  bad  but  the  differences 
from  a  multi-layered  model  are  very  interesting.  The  comparisons  of  the  mid-  depth 
currents  off  Somalia  with  observations  is  very  rewarding.  Once  again  I  hope  that 
this  work  convinces  you  that,  if  we  use  ocean  models  with  a  shape  like  the  actual 
ocean,  a  forcing  estimated  from  real  atmospheric  data  and  very  fine  grid  spacing, 
we  get  marvellous  results  that  convince  us  that  we  are  learning  about  the  oceans 
from  our  calculations. 


James  J.  O’Brien 
Director 

Mesoscale  Air-Sea  Interaction  Group 


Abstract 


A  new  numerical  ocean  model  with  multiple  isopycnal  layers,  has  been  used 
to  model  the  Indian  Ocean.  Normal  vertical  modes  are  used  for  initialization  and 
in  a  new  open  boundary  formulation.  A  21  year  integration  with  the  Hellerman- 
Rosenstein  wind  stress  is  made  with  a  3.5  layer  and  a  1.5  layer  version  of  the  model. 

The  solution  with  three  active  layers  reproduces  the  observed  general  cir¬ 
culation  and  variability  of  the  Indian  Ocean,  for  instance^the  semi-annual  equatorial 

undercurrent  and  Yanai  wave  field  in  the  west.  The  seasonal  changes  in  the  Somali 
» 

Current  system  /s  studied  in  more  detail.  It  is  found  that  barotropic  instability  is 
likely  to  cause  the  generation  of  the  Great  Whirl  in  early  June.  We  find  a  very 
good  agreement  between  the  observed  undercurrents  and  the  simulations  in  the 
model.  Equatorial  onshore  flow  below  the  thermocline  in  June  is  associated  with 
the  disappearence  of  the  undercurrent  below  the  Somali  current.  The  return  of 
this  undercurrent  in  the  fall  is  caused  by  instability  of  the  Great  Whirl.  Experi¬ 
ments  where  the  duration  of  the  summer  monsoon  is  extended  show  that  the  initial 
decrease  in  the  magnitude  of  the  Great  Whirl  is  due  to  eastward  and  downward 
energy  transfer  rather  that  due  to  relaxation  of  the  wind.  The  solutions  of  the 
model  indicate  that  baroclinic  instability  plays  an  important  role  in  the  decay  of 
the  Great  Whirl. 

The  solution  with  a  single  active  layer  is  essentially  the  same  for  the  upper 
layer  until  in  the  late  summer  monsoon,  when  the  flow  becomes  unstable.  Different 
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decay  patterns  of  the  whirl  and  associated  eddies  leads  to  different  flows  during  the 
winter  monsoon. 
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1.  INTRODUCTION 

1.1.  Motivation  and  objectives 

There  has  been  a  long  tradition  at  Florida  State  University  of  using  isopy- 
cnic  layered  models  with  a  limited  vertical  resolution  to  study  ocean  dynamics.  Mc- 
Nider  and  O’Brien  (1973)  used  a  four-layer  model  to  study  coastal  upwelling,  while 
Hurlburt  and  Thompson  (1980)  investigated  the  circulation  in  the  Gulf  of  Mex¬ 
ico  using  three  different  models:  a  two-layer,  a  reduced-gravity,  and  a  barotropic 
model.  Studies  using  reduced-gravity  models  with  one  moving  layer  have  been  suc- 
cesful  in  simulating  equatorial  flow,  for  instance  El  Nino  in  the  Pacific  (Busalacchi 
and  O’Brien,  1980).  The  work  by  Hurlburt  and  Lin,  (1981),  who  demonstrated 
that  a  reduced-gravity  model  contained  the  essential  physics  for  generation  of  a  So¬ 
mali  Current,  motivated  later  work  with  realistic  geometries  for  the  Indian  Ocean 
as  mentioned  below,  e.g.,  Luther  and  O’Brien  (1985).  These  studies  suggests  that 
further  progress  can  be  made  by  including  additional  physics.  For  instance,  Zebiak 
and  Cane  (1987)  added  a  mixed  layer  of  constant  depth  with  thermodynamics  on 
top  of  the  upper  layer  and  coupled  the  ocean  to  an  atmospheric  layer  to  study  the 
El  Nino-Southern  Oscillation  phenomenon.  In  a  similar  study  (Schopf  and  Suarez, 
1988)  a  two-layer  reduced-gravity  model  allowed  temperature  variations  in  both 
active  layers. 

Layers  permit  accurate  representation  of  low  vertical  modes  compared  to 
fixed  levels.  Therefore  for  fixed  resources  one  can  have  finer  horizontal  resolution. 
Bleck  and  Boudra  (1981,  1986)  compared  results  from  quasi-isopycnic  layered  mod- 
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els  with  other  numerical  formulations,  for  instance  quasi-geostrophic  layered  models 
and  isobaric  models  as  Bryan’s  (1969)  model.  They  found  that  the  latter  two  sup¬ 
pressed  barotropic  and  baroclinic  instability.  Two  and  three  layer  versions  of  their 
quasi-isopycnic  model  were  also  used  to  study  the  effect  of  model  parameters  on 
the  circulation  of  the  South  Atlantic-Indian  Ocean  (Boudra  and  Ruijter,  1986). 

The  principal  advantage  of  using  density  as  the  vertical  coordinate  as  in 
the  models  above,  compared  to  depth,  is  that  no  artificial  cross  isopvcnal  mixing 
occurs.  Another  consequence  of  using  the  layer  thicknesses  as  variables  is  that  bet¬ 
ter  vertical  resolution  automatically  is  obtained  in  areas  with  strong  stratification 
(Bleck  and  Boudra,  1986).  The  low  vertical  resolution  effectively  filters  out  higher 
vertical  modes,  which,  due  to  their  smaller  length  scales,  are  not  adequately  resolved 
by  the  horizontal  discretization  in  the  numerical  model  and  consequently  produce 
unwanted  noise.  Finally,  these  models  with  coarse  vertical  resolution  and  without 
prognostic  thermohaline  equations  require  relatively  small  computational  resources 
compared  to  oceanic  general  circulation  models  with  full  physics.  This  makes  it 
possible  to  perform  numerical  experiments  where  the  effect  of  changing  various 
parameterizations,  boundary  conditions,  geometries  and  wind  stress  fields  can  be 
investigated.  The  multi-layer  formulation  to  be  used  in  this  study  is  described  in 
section  2. 

The  purpose  of  this  study  is  to  build  a  new  multi-layer  ocean  model,  which 
can  simulate  the  seasonal  changes  of  surface  and  subsurface  currents  in  the  North¬ 
west  Indian  ocean,  in  particular  the  Somali  Current  system.  As  forcing,  a  seasonal 
climatological  wind  stress  is  used  to  obtain  a  quasi-periodic  ocean  circulation.  The 
model  results  will  be  analyzed  and  compared  with  observations.  The  effect  of  addi¬ 
tional  layers  will  also  be  studied  by  comparing  the  currents  in  the  upper  layer  of  the 
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new  model  with  calculations  from  a  reduced-gravity  model  with  the  same  geometry 
and  forcing. 

1.2.  The  Somali  Current  System 

The  circulation  of  the  northwestern  part  of  the  Indian  Ocean  is  highly  in¬ 
fluenced  by  the  monsoon  winds.  The  Somali  Current,  which  has  a  volume  transport 
comparable  to  that  of  the  Gulf  Stream  (Lighthill,  1969),  changes  direction  with  the 
Monsoon  winds:  Northeastward  currents  during  the  summer  (Southwest  Monsoon) 
and  southwestward  currents  during  the  winter  (Northeast  Monsoon).  Swallow  and 
Bruce  (1966)  reported  a  transport  of  60  Sv  in  the  upper  200  m  during  the  summer 
monsoon,  and  surface  velocities  up  to  3.7  m/s  have  been  measured  (Duing  et  ah, 
1980).  Measurements  made  in  1964  during  the  International  Indian  Ocean  Expe¬ 
dition  (IIOE)  and  during  the  seventies  showed  that  the  flow  pattern  of  the  Somali 
Current  is  far  more  complicated  than  the  continuous  northward  or  southward  flow 
as  seen  in  climatological  atlases,  e.g.  Schott  (1983).  In  particular,  during  the  sum¬ 
mer  it  is  now  anticipated  that  a  two-gyre  system  exists.  Knox  and  Anderson  (1985) 
defines  the  Somali  Current  System  as  “the  currents  along  the  African  coast  plus 
their  associated  eddies,  offshore  meanders  and  recirculations.” 

1.3.  Observations 

The  observed  seasonal  changes  have  been  reviewed  by  Schott  (1983),  Knox 
and  Anderson  (1985)  and  can  be  summarized  as  follows:  In  late  winter  the  Somali 
current  has  a  deep  southward  flow  from  the  island  of  Socotra  to  about  2°-3°  S. 
Here  the  flow  joins  the  East  African  Coast  current  coming  from  the  south  into 
the  eastward  Equatorial  Counter  Current.  In  the  early  spring,  usually  in  March,  a 
northward  surface  current  is  established  north  of  5°N  due  to  local  wind  stress  curl. 
Below,  between  150  m  -  600  m,  a  southward  undercurrent  is  found  (Quadfasel 
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and  Schott,  1983).  Between  the  equator  and  5°N  the  surface  flow  is  southward. 
South  of  the  equator  a  northward  current  increases  in  strength  at  this  time,  crosses 
the  equator  and  flows  offshore.  During  April  and  May  this  current  migrates  to 
the  north.  Southwesterly  winds  start  in  early  May,  and  strong  upwelling  is  found 
north  of  3°-5°  N,  which  is  the  latitude  where  offshore  flow  is  found.  A  southward 
return  flow  develops  offshore,  and  a  gyre,  the  so  called  Southern  Gyre,  is  formed. 
As  the  coast-parallel  wind  strengthens  and  strong  anticyclonic  wind  stress  curl 
occurs  offshore,  a  northern  gyre,  often  called  the  Great  Whirl,  a  name  introduced 
by  Findlay  (1866),  who  first  reported  it,  is  formed  between  5°  and  10°  N.  Wedges 
of  cold  water  are  found  along  the  coast  north  of  each  gyre.  As  the  Great  Whirl 
deepens,  the  southward  undercurrents  disappear.  This  two  gyre  system  is  stable 
until  August  or  September,  when  the  Southern  Gyre  propagates  northward  and 
merges  with  the  Great  Whirl,  as  observed  by  Bruce  (1973).  The  salinity  of  the 
Southern  Gyre  is  lower  than  in  the  surrounding  water  mass  before  the  coalescense. 
suggesting  mass  transport  from  the  south.  After  this  event  the  Somali  Current 
becomes  stationary,  flowing  from  4°  S  to  10°  N,  with  northward  transport  above 
150  m  and  southward  transport  reappears  below  (Quadfasel  and  Schott,  1983).  To 
the  north-east  of  the  Great  Whirl  a  third  warm  core  eddy,  the  Socotra  eddy  is  found 
in  the  late  summer.  Its  relative  high  salinity  indicates  that  the  water  in  this  eddy 
is  advected  from  a  higher  latitude.  This  situation  lasts  until  after  the  onset  of  the 
northeast  monsoon  in  November,  when  a  shift  to  southward  currents  takes  places, 
first  along  the  coast  and  later  offshore.  During  the  winter  monsoon  the  flow  is  a 
deep  southward  continuous  flow  from  Socotra  to  5°S. 

The  seasonal  cycle  as  described  above  is  considered  the  most  common 
situation  and  was  the  situation  during  the  FGGE  year  1979.  Swallow  and  Fieux, 
(1982)  reviewed  historical  data  for  May  and  June,  from  1900  to  1973.  They  found 
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that  the  two  gyre  situation  was  the  most  common.  The  southern  gyre  was  missing 
in  three  years,  and  could  be  found  during  May  or  June  in  30  of  the  years.  Out  of 
69  years,  evidence  of  the  Great  Whirl  was  found  in  55. 

1.4.  Modelling 

It  was  suggested  by  Lighthill  (1969)  that  the  onset  of  the  Somali  Current 
in  the  spring  was  a  result  of  remote  forcing  due  to  winds  over  the  Arabian  Sea 
turning  eastwards.  The  linear  theory  predicted  a  response  time  of  one  month  for 
the  coastal  current  to  be  set  up  associated  with  the  reflection  of  westward  travelling 
baroclinic  equatorial  Rossby  waves.  The  theory  also  predicted  the  existence  of  a 
deep  undercurrent,  and  an  offshore  countercurrent,  where  the  magnitude  of  the 
undercurrent  depends  on  the  relative  strength  of  the  first  baroclinic  and  barotropic 
modes.  His  most  important  result  was  the  short  time  scale  required  to  create  a 
western  boundary  current  near  the  equator  compared  to  the  long  time  needed  at 
mid-latitudes. 

However,  Leetmaa  (1972,  1973)  observed  that  in  1970  and  1971,  a 
northerly  flow  forced  by  local  winds  started  in  the  south,  off  the  coast  of  Kenya, 
some  time  before  the  onset  of  the  westerlies  in  the  interior.  This  raised  a  contro- 
versary  whether  the  Somali  Current  was  a  result  of  local  or  remote  forcing. 

Lin  and  Hurlburt  (1981)  used  a  simple  reduced-gravity  model  in  a  rectan¬ 
gular  domain  to  model  the  response  to  local  meridional  wind  forcing.  They  found  a 
response  time  of  one  week,  and  that  a  sequence  of  cyclonic  and  anticyclonic  eddies 
formed  north  of  the  equator  and  moved  poleward.  Their  simple  model  produced 
the  most  important  features:  Coastal  upwelling  was  seen  to  the  northwest  of  the 
latitude  where  the  currents  of  a  large  warm  eddy  turned  offshore,  and  found  and 
wedges  of  a  thin  upper  layer  500  km  offshore  just  north  of  that  separation  point. 
The  more  detailed  solution,  i.e.,  the  magnitude  of  the  transport,  eddy  activity  and 
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intensity  of  the  upwelling  depended  on  whether  open  or  closed  boundary  conditions 
were  used  and  of  the  width  of  the  wind  stress  field. 

The  first  numerical  simulation  of  the  Indian  Ocean  circulation  is  due  to 
Cox  (1970).  The  model  had  realistic  coastlines,  with  the  horizontal  resolution  of 
1°  and  7  layers  in  the  vertical.  The  observed  large  scale  structure  of  the  Indian 
Ocean  was  reproduced,  including  the  annual  reversal  of  the  Somali  Current  with 
coastal  upwelling  during  the  summer,  but  eddies  like  the  Great  Whirl  and  the 
southern  gyre  were  missing  due  to  the  coarse  horizontal  resolution  in  the  model. 
During  the  summer  the  model  boundary  current  did  not  extend  as  far  north  as 
observed,  and  the  velocities  in  the  boundary  current  were  a  factor  of  two  or  more 
too  small.  It  was  suggested  that  the  Somali  Current  is  barotropic  up  to  3°N,  where 
coastal  upwelling  drive  a  baroclinic  transport.  Using  Ekman  theory,  Cox  (1970) 
found  vertical  velocities  up  to  5  m/day  due  to  local  winds.  Since  these  winds  are 
important  for  the  upwelling,  it  was  concluded  that  local  forcing  was  important  in 
the  north,  while  the  southern  part  is  a  remotely  forced  western  boundary  current 
in  a  classical  sense,  he.,  a  return  flow  in  response  to  a  Sverdrup  interior. 

In  a  later  study,  Cox  (1976)  examined  the  two  theories  of  remote  versus 
local  forcing  using  the  numerical  model  of  Cox  (1970),  but  applied  it  to  a  large 
rectangular  ocean,  symmetric  around  the  equator.  It  was  found  that  the  local 
forcing  was  dominant  initially  and  caused  the  current  to  extent  further  north  than 
by  the  remote  forcing,  which  became  important  after  a  couple  of  weeks.  A  linear 
analytical  tneory  by  Anderson  and  Rowlands  (1976)  corroborates  these  results, 
showing  that  the  amplitude  of  the  boundary  current  increases  linearly  in  time  in 
the  cause  of  local  wind  forcing,  but  quadratic  in  time  for  remote  forcing. 

Using  monthly  means  of  climatological  wind  stress  to  force  a  reduced- 
gravity  model,  Luther  and  O’Brien  (1985)  and  Luther  et  al.,  (1985)  succesfully 
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reproduced  most  of  the  observed  features  of  the  seasonal  cycle  as  described  above. 
Their  simulation  was  the  first  eddy  resolving,  fully  nonlinear  model  with  realistic 
geometry  of  the  northwestern  Indian  Ocean.  During  and  after  the  collapse  of  the 
two  gyre  system  some  differences  between  model  and  observations  are  seen,  and 
of  course  their  model  does  not  contain  information  of  the  vertical  structure  of  the 
currents. 

The  scenario  described  in  the  section  above,  is  subject  to  some  interannual 
variability.  Simmons  et  al.,  (1988)  found  that  the  number  of  eddies  around  Socotra 
in  the  fall  depended  on  the  wind  forcing,  which  was  based  on  ship  observations  from 
different  years.  In  a  recent  study,  Luther  and  O’Brien  (1989),  also  using  winds  from 
ship  observations  applied  to  the  model  mentioned  above,  found  that  out  of  23  years 
of  model  simulations,  the  southern  gyre  was  missing  in  two.  The  coalescense  of  the 
southern  gyre  and  the  Great  Whirl  was  seen  in  all  but  7  years.  On  the  other  hand, 
the  Great  Whirl  and  the  Socotra  eddy  occurred  every  summer. 

Modelling  the  undercurrents  has  been  less  succesful.  Hurlburt  and  Thomp¬ 
son  (1976)  found  a  cyclonic  inflow  and  an  undercurrent  associated  with  the  gen¬ 
eration  of  the  Great  Whirl  in  a  two-layer  model.  The  upper  layer  response  was 
in  good  agreement  with  the  observed  whirl  and  its  northward  migration,  but  the 
flow  in  the  deep  layer  was  not.  Using  a  linear  vertical  mode  model  McCreary  and 
Kundu  (1985)  found  that  an  undercurrent  developed  south  of  an  area  with  wind 
stress  curl.  For  winds  without  curl,  no  undercurrent  was  seen.  In  a  later  study, 
McCreary  and  Kundu  (1988)  used  a  two-layer  reduced-gravity  model  and  found 
that  no  undercurrent  was  produced  even  in  the  presence  of  wind  stress  curl.  Schott 
(1987)  compared  observations  with  the  output  from  Philanders’s  27-level  model 
of  the  Indian  Ocean.  The  results  are  not  published,  but  the  model  is  described 
in  Philander  and  Pacanowski  (1984,  1986).  The  model  was  driven  by  winds  from 
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Hellerman  and  Rosenstein  (1983)  and  initiated  with  climatological  temperature  and 
salinity  data  given  by  Levitus  (1982).  A  good  qualitative  agreement  between  ob¬ 
servations  and  model  was  found  for  the  Somali  Current.  The  model  reproduced  the 
undercurrent  during  the  North  East  Monsoon,  but  lacked  a  deep  response  found  in 
observations.  The  latter  might  be  due  to  the  rather  short  spin  up  time  of  3  years 
from  initialization  of  the  model.  Other  recent  modelling  efforts  are  summarized  by 
Luther  ( 1987). 

1.5.  Model  Experiments 

Schott  (1987)  lists  several  topics  not  yet  fully  investigated  due  to  lack  of 
observations  and  model  results: 

Most  emphasis  has  been  given  to  the  summer  monsoon,  and  little  is  known 
about  the  structure  of  the  currents  during  the  winter  monsoon.  The  extent  of  the 
undercurrents  is  not  yet  known.  Other  questions  to  be  addressed  are  whether 
barotropic-baroclinic  instability  is  important  for  the  collapse  of  the  Great  Whirl, 
and  whether  local  wind  changes  may  be  responsible  responsible.  The  cause  of 
this  break  down  is  not  yet  understood.  These  problems  will  be  addressed  here,  in 
particular  the  two  latter  topics. 

Two  different  configurations  of  a  new  numerical  model  will  be  used  in  this 
study.  A  version  with  one  active  layer  over  a  non-moving  lower  layer  will  produce 
results  which  has  the  same  physics  as  in  the  models  by  Luther  and  O’Brien  (1985), 
Luther  et  al.  (1985),  but  is  driven  by  a  different  wind  stress.  Here  a  pseudo  wind 
stress  based  on  the  Hellerman  and  Rosenstein  (1983)  wind  stress  analysis  is  used  as 
wind  forcing.  The  investigations  mentioned  above  used  another  climatological  wind 
stress,  e.g.,  the  NOAA  Global  Marine  Sums  data  set.  The  second  model  version  will 
have  three  moving  layers  over  an  infinitely  deep  lower  layer.  The  geometry,  wind 
forcing,  and  physical  parameters,  with  the  exception  of  vertical  density  profile,  is 
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not  changed  from  the  computations  with  a  single  layer,  in  order  that  the  effect  of 
additional  vertical  resolution  can  be  studied. 

The  models  are  spun  up  for  21  years,  and  an  analysis  of  their  seasonal 
cycle  is  presented  and  compared  with  observations.  First  the  features  of  the  general 
circulation  of  the  Indian  Ocean  are  described  using  10  year  averages  of  the  model 
results  from  year  1 1  to  20.  Next  we  focus  on  the  seasonal  variability  of  these  averages 
for  the  Somali  Current  system  including  the  undercurrents.  Finally  the  evolution 
during  the  summer  monsoon  in  year  21  is  described  and  analysed  in  detail  in  order 
to  determine  the  dynamics  involved  in  the  decay  of  the  Great  Whirl. 


2.  MODEL  FORMULATION 


2.1.  Physics  of  the  model 

We  apply  the  equations  of  motion  and  conservation  of  mass  for  the  ocean 
in  spherical  coordinates  (see  Semtner,  1986).  Let  the  longitude  and  latitude  be 
given  by  <f>  and  $  and  the  velocity  components  towards  the  east  and  north  be  u 
and  v,  respectively.  As  radial  coordinate  we  use  z  =  r  —  a,  where  r  is  the  radial 
distance  from  the  center  of  a  spherical  earth  with  radius  a.  We  choose  2  =  0  to 
be  the  surface  of  the  ocean  at  rest.  Define  vertically  integrated  volume  transport 
components  U j  and  Vj  by 
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between  two  surfaces  6,t)  and  Zj+i(<t>,  0,  t),  with  an  equivalent  expression  for 
Vj.  The  thickness  of  the  j-th  layer  defined  by  this  integration,  is  Hj  =  (zj+\  —  Zj). 

Consider  an  ocean  consisting  of  several  layers  of  uniform  density  as  shown 
in  Fig.  1.  The  layers  are  labelled  with  increasing  numbers  downward.  Let  us  assume 
that  all  layers  have  a  positive  thickness  everywhere  for  all  time.  This  implies  that 
layers  are  not  allowed  to  surface  or  merge,  and  that  the  bottom  topography,  given 
by  z  =  D(<f>,0),  is  always  in  the  lowest  layer.  A  numerical  technique,  the  Flux 
Corrected  Transport  (FCT)  scheme  (see  Zalesak,  1979),  makes  it  possible  to  relax 
these  restrictions  and  has  recently  been  applied  to  layered  ocean  models  to  allow 
fronts  at  the  surface,  e.g.,  Bleck  and  Boudra  (1986),  Huang  (1987).  The  method 
requires  that  the  flux  in  the  continuity  equation  is  calculated  as  a  combination 
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I 

^  Figure  1.  Vertical  structure  of  a  four  layer  isopycnal  model.  Bottom  topography 

does  not  intersect  isopycnals. 
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from  a  diffusive  scheme  and  a  non-diffusive  scheme  with  the  constraint  that  the 
layer  thickness  always  is  greater  than  zero.  One  problem  with  the  method  is  that 
it  does  not  allow  quadratic-conservative  flux  form  of  the  non-linear  terms,  which 
requires  division  by  the  layer  thickness.  Consequently  it  will  not  be  applied  in  the 
present  study.  The  derivation  of  the  equations  below  is  given  in  Appendix  A.  If  the 
density  and  thickness  of  the  j-th  layer  is  given  by  pj  and  Hj,  the  transport  equation 
for  U  becomes 


al±  +  1  a  (vl\  .  \  i  (W\  _  IMi  _  fV. 

dt  acosO  d4>\Hj)  a  dO  \  Hj  /  aHjCotO  3 


-Hj 
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dt 
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dr,  U,  JlH, 

+  Wgi  -  9  2^(Pj  ~  Pi) 


l—l 


d(f>  J 
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where  the  horizontal  friction  for  the  U  equation  is  given  by 
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3  \HjJ  a2  cos20 


Uj(l  -2cos20)  +  2sin0Hj^(^)] 


and  similarly  for  V,  we  have 
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with  the  horizontal  friction  term 
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(2-5) 


In  the  equations  above  g  is  the  acceleration  of  gravity,  r  the  tangential 
stress  due  to  vertical  friction,  where  the  superscripts  denote  the  <f>  or  0  component 
and  top  or  bottom  of  layer.  With  N  layers  the  surface  displacement  rj  is  given  by 

N 

T,  =  D  +  jTH, 

1=1 


(2-6) 
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In  the  equation  above  the  parameterization  of  the  horizontal  eddy  viscosity  is  based 
on  the  velocity  field  and  not  the  vertically  integrated  transport  as  most  often  seen  in 
layered  models.  The  formulation  has  been  chosen  so  a  thin  layer  with  large  velocities 
will  experience  more  friction  than  a  thick  layer  with  equivalent  transport.  It  is  also 
seen,  that  if  the  layer  thickness  gradient  is  constant,  the  geostrophic  current  will 
experience  a  spin  down  due  to  friction.  In  the  case  where  the  horizontal  friction 
is  proportional  to  the  Laplacian  of  the  transport,  a  steady  geostrophic  solution 
is  possible.  The  magnitude  of  the  horizontal  friction  coefficient,  A ,  used  in  the 
calculations  was  750  m^/s,  the  same  value  used  by  Woodberry  et  al.  (1989). 

The  only  vertical  stress  to  be  applied  in  the  model  in  this  study  is  the 
wind  stress,  which  acts  as  a  body  force  on  the  upper  layer. 


The  continuity  equation  becomes 

OH;  1  r  dU<  d  1 

~W  +  +  dol'’’cos0}}  = 


(2-7) 


where  we  is  a  source  term  due  to  entrainment.  This  term  is  positive  for  the  upper 
layer  in  case  its  thickness  becomes  less  that  a  preset  minimum  depth  Hmjn,  and  of 
the  same  magnitude  but  negative  for  the  second  layer.  For  deeper  layers  the  term 
is  always  set  to  zero.  This  entrainment  is  included  only  to  prevent  the  interface 
between  the  first  two  layers  to  surface.  The  parameterization  of  McCreary  and 
Kundu,  (1988)  is  used.  However,  the  effect  on  the  upper  layer  density,  momentum 
and  kinetic  energy  balance  by  this  entrainment  has  been  neglected  here.  For  the 
momentum  equations  the  term  ignored  is  wev 2,  which  is  the  transfer  of  momentum 


from  the  second  to  the  top  layer.  The  entrainment  velocity  is  given  by 


(2-8) 


(Hi-Hnin) 


1\IJ  1  — 11  min )  it  ^  IT 

-  Umin  (2  _  8) 

0  H\  >  Hmin 

The  simplification  made  by  only  considering  the  effect  of  we  in  the  conti¬ 


nuity  equation  corresponds  to  entrainment  where  the  engulfed  water  into  the  upper 
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layer  has  zero  velocity  and  is  heated  instantaneously  to  the  density  p\.  It  should 
be  emphasized  that  this  adjustment  is  rarely  active  in  the  results  presented  in  this 
work.  We  chose  60  m  for  Hmin  and  a  time  constant  re  of  1.2  hours. 

We  can  ignore  forcing  due  to  the  gradient  in  atmospheric  pressure,  pa,  com¬ 
pared  to  the  wind  stress  for  large  scale  and  mesoscale  motion,  e.g.  Gates  (1966); 
thus  the  barotropic  pressure  gradient  is  contained  in  the  term  where  the  surface 
deviation  tj  appears.  Because  of  the  large  phase  speeds  of  barotropic  gravity  waves 
most  numerical  models  have  a  special  treatment  of  this  mode.  For  instance,  most 
authors  filter  out  these  waves  by  applying  the  rigid  lid  approximation,  for  example 
as  in  Bryan’s  (1969)  world  ocean  model,  or  using  a  semi-implicit  numerical  scheme 
(O’Brien  and  Hurlburt,  1972,  Hurlburt,  1974,  Hurlburt  and  Thompson,  1976).  Here 
we  shall  apply  another  method  to  remove  the  barotropic  modes,  including  the  plan¬ 
etary  waves,  by  assuming  that  the  pressure  gradient  vanishes  in  the  lowest  layer, 
which  implies  that  the  velocity  also  is  zero  in  the  deep  ocean.  From  (2-2)  and  (2-4) 
we  obtain  by  taking  the  limit  H n  — »  oo,  that  the  gradient  of  the  surface  elevation 
is  given  by 

Vri=  V  (PN  ~~-Pi)vHi  (2-9) 

i^\KpN 

The  effect  of  bottom  topography  can  to  some  extent  be  included  when  (2  —  9)  is 
applied  by  generalizing  the  method  of  Cushman-Roisin  and  O’Brien  (1983).  They 
demonstrated  that  the  effect  of  variable  bottom  topography  in  a  two-layer  model  can 
be  simulated  by  locally  changing  the  phase  speed  in  a  reduced-gravity  model.  Here 
we  will  use  (2  —  9)  with  N  =  4  and  N  =  2,  which  gives  us  a  3.5  layer  model  and  a  1.5 
layer  model,  respectively.  Explicit  integration  of  the  set  of  equations  (2  —  2)  to  (2  —  8) 
with  a  finite  depth,  providing  that  the  isopycnals  does  not  intersect  the  bottom 
topography  (e.g.,  Fig.  1),  is  in  principle  done  in  the  same  way,  but  limitations  due 
to  the  very  short  time  step  in  the  integration  makes  this  impractical.  The  computer 
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time  needed  increases  by  a  factor  of  40  compared  to  the  reduced-gravity  case  with 
the  same  number  of  layers. 

2.2.  Model  geometry 

The  model  covers  the  northwestern  part  of  the  Indian  Ocean  from  25.1°S 
to  26.1°N  and  34.8°E  to  119.6°E.  The  coastlines  are  identical  to  those  used  by 
Woodberry.  et  al..  (1989).  As  seen  from  (Fig.  2)  the  eastern  and  southern  boundaries 
are  open,  while  the  no-slip  condition  is  applied  along  land  boundaries.  The  200  m 
isobath  was  used  to  define  the  coastlines.  This  implies  that  groups  of  coral  reef 
islands,  such  as  the  Laccadives  and  Maldives  southwest  of  India,  and  the  Seychelles 
Bank,  Sava  de  Malha  Bank  with  the  Agalegas  Islands  and  Nazarath  Bank  with 
the  Cargados-Carajos  Islands,  which  am  fjund  northeast  and  west  of  Madagascar, 
appears  as  large  islands. 

It  might  have  been  preferable  to  extend  the  model  further  south  to  include 
more  of  the  southern  gyre,  and  model  the  flow  around  the  southern  tip  of  Madagas¬ 
car,  but  the  present  geometry  was  chosen  since  prepared  wind  fields  and  digitized 
coastlines  were  available. 

2.3.  Numerical  formulation 

The  model  equations  are  discretized  in  space  on  the  Arakawa  C-grid 
(Mesinger  and  Arakawa,  1976,  Arakawa  and  Lamb,  1977).  The  numerical  scheme, 
which  is  mass  and  energy  conserving,  is  nearly  identical  to  the  C  scheme  tested  by 
Grammeltvedt  (1969)  and  first  proposed  by  Lilly  (1965).  With  a  distance  between 
two  similar  points,  for  instance  two  //-points,  of  0.2°  in  both  horizontal  directions, 
the  model  domain  contains  425  x  256  grid  points  per  layer  for  each  variable. 

With  the  barotropic  mode  excluded  using  (2-9)  we  can  solve  the  finite 
difference  equivalents  of  eqns.  (2-2)  to  (2-8)  using  an  explicit  time  integration  scheme 
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INDIAN  OCERN  MODEL 


ERST  LONGITUDE 


Figure  2.  Model  geometry.  The  200  m  isobath  is  used  as  coastlines.  There  are 
open  boundaries  to  the  south  and  in  the  east. 
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without  a  severe  restriction  in  time  step.  With  the  grid  above  the  leap-frog  time 
integration  scheme  is  stable  with  a  time  step  of  about  20  minutes,  depending  on 
the  current  speed  and  the  phase  speed  of  the  first  vertical  mode  in  the  model.  For 
the  Laplacian  friction  term  a  Dufort-Frankel  implicit  scheme  is  applied  for  stability 
(e.g.,  O'Brien,  1986).  An  Euler  forward  scheme  is  applied  to  every  99  steps  to  filter 
out  the  computational  mode  inherent  in  the  leap-frog  scheme. 

The  model  is  coded  in  Fortran  200,  fully  vectorized,  to  run  on  the  Cyber 
205  and  ETA  10  vector  computers.  Less  than  5%  of  the  CPU  time  is  spent  for  scalar 
computations.  Half  precision  (32  bits)  variables  have  been  used  in  the  calculations, 
to  make  the  code  nearly  twice  as  fast  compared  to  when  full  precision  (64  bits) 
variables  are  declared.  The  largest  terms,  i.e.,  Coriolis  terms,  pressure  gradient 
terms  and  wind  stress  are  computed  separately  to  avoid  unnessessary  round  off 
errors.  Their  sum  is  added  to  the  sum  of  the  remaining  terms  in  the  momentum 
equation.  For  the  same  reason,  only  the  layer  thickness  deviations,  (Hj  —  Hqj ),  was 
stored  and  used  as  dependent  variable.  The  model  requires  about  25  minutes  of 
CPU  time  per  active  layer  per  model  year  on  a  one  processor  ETA10-G  or  1  hour 
on  the  Cyber  205.  The  memory  requirement  for  the  ocean  model  is  designed  to  fit 
within  a  4  megaword  central  memory  to  avoid  paging.  The  3.5  layer  model  requires 
2.4  megawords  of  storage. 

2.4.  Boundary  and  initial  conditions 

At  the  coast  it  is  assumed  that  both  components  of  the  transport  vanishes. 
i.e.,  the  no  slip  condition.  The  boundary  layer  is  not  resolved,  but  as  shown  by  Cox 
( 1979),  this  condition  is  provides  a  vorticity  flux  from  the  boundary  into  the  interior, 
which  gives  a  more  realistic  solution.  While  these  boundary  conditions  are  simple, 
the  correct  conditions  to  be  used  at  open  boundaries  in  multi-layer  models  are  less 
obvious.  An  easy  approach  is  to  apply  a  Sommerfeld  radiation  condition  to  each 
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layer  separately.  However,  different  vertical  modes  have  different  phase  speeds,  and 
a  more  accurate  method  is  to  separate  the  solution  along  the  open  boundary  into 
these  modes.  Using  linear  theory  for  an  ocean  with  flat  bottom,  Lighthill  (1969) 
showed  that  the  vertical  modes  of  the  velocity  are  eigenvectors  to  the  matrix 

v  =  *=?**!  Hk  (2  -  10) 

Pj 

and  the  eigenvalues  correspond  to  the  equivalent  depths  for  each  mode.  By  inverting 
the  matrix  which  has  the  eigenvectors  of  aj j.  as  columns,  we  can  find  the  amplitude 
for  each  vertical  mode  from  the  solution  of  the  layered  model.  The  radiation  condi¬ 
tion  is  applied  to  each  mode  separately  and  the  resulting  current  components  along 
the  open  boundary  can  then  be  computed  as  the  sum  of  these  vertical  modes. 

Several  numerical  approximations  to  the  Sommerfeld  radiation  condition 
have  been  proposed.  The  scheme  implemented  in  this  model  was  given  by  Camer- 
lengo  and  O’Brien  (1980),  who  improved  and  simplified  the  original  scheme  by 
Orlanski  (1976).  Both  these  conditions  and  the  more  simple  extrapolation  condi¬ 
tion  were  tested  on  a  model  problem  to  select  the  most  suitable  open  boundary 
condition,  and  the  best  results  were  obtained  using  the  Camerlengo  and  O’Brien 
scheme  on  each  vertical  mode.  Details  are  given  in  Appendix  B.  A  review  of  radia¬ 
tion  boundary  conditions  can  be  found  in  Hedley  and  Yau  (1988),  while  R.0ed  and 
Cooper  (1986,  1987)  considered  a  wider  range  of  open  boundary  conditions.  In  the 
three  references  above  it  was  demonstrated  that  the  schemes  based  on  the  radiation 
condition  may  be  inadequate  in  some  cases. 

The  vertical  modes  for  the  Indian  Ocean  given  by  Gent,  et  al.,  (1983) 
have  been  used  to  select  the  initial  layer  thicknesses.  For  the  first  three  vertical 
baroclinic  modes  they  reported  equivalent  depths  of  79.9  cm,  30.5  cm  and  12.6  cm, 
respectively.  By  selecting  pj  for  four  layers  as  1.0239,  1.0262,  1.0273  and  1.0279 
g/cm'5  and  layer  thicknesses  II j,  200  m,  250  m  and  400  m  for  the  three  upper  layers, 
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we  obtain  equivalent  depths  of  102  cm,  22.5  cm  and  10.2  cm.  This  corresponds  to 
Kelvin  wave  speeds  of  316  cm/s,  149  cm/s  and  100  cm/s  for  the  three  baroclinic 
modes,  compared  to  280  cm/s,  173  cm/s  and  111  cm/s  found  by  Gent,  et  ah,  (1983). 
Realistic  initial  phase  velocities  for  the  internal  gravity  waves  were  chosen  to  ensure 
good  phase  correlation  with  observed  currents.  For  comparison,  the  Kelvin  wave 
phase  speeds  of  the  first  three  vertical  modes  in  the  model  by  Cox  (1976)  were  291 
cm/s,  179  cm/s  and  117  cm/s.  In  the  2.5  layer  model  of  McCreary  and  Kundu 
(1988)  the  initial  values  were  321  cm/s  and  123  cm/s.  For  the  1.5  layer  model 
pi  =  1.025  g/cm'5  and  pi  =  1.028  g/cm^  was  selected,  which  corresponds  to  a  value 
of  reduced  gravity  of  0.03  m/s2  also  used  by  Woodberry,  et  al.,  (1989).  This  results 
in  a  Kelvin  wrave  speed  of  245  cm/s.  Kindle  and  Thompson  (1989)  used  the  same 
density  difference,  but  used  an  initial  upper  layer  thickness  of  250  m  in  their  1.5 
layer  model,  resulting  in  a  phase  speed  of  274  cm/s. 

However,  when  the  models  are  fully  non-linear  as  those  presented  here, 
the  actual  phase  speed  will  vary  with  space  and  time  depending  of  the  solution. 

2.5.  Wind  forcing 

The  model  is  forced  by  a  climatological  monthly  mean  wind  stress  based 
on  the  data  set  prepared  by  Hellerman  and  Rosenstein  (1983).  A  pseudo-stress  is 
formed  by  dividing  these  data  by  the  product  of  an  average  drag  coefficient  and 
an  air  density  (Woodberry,  et  ah,  1989).  This  allows  the  drag  coefficient  and  air 
density  to  be  model  parameters  independent  of  the  wind  analysis.  A  constant  drag 
coefficient  of  1.5  10~^  and  an  air  density  of  1.2  kg/m^  has  been  applied  here.  Figs.  3 
to  6  show  the  wind  stress  and  its  curl  for  February,  May,  August  and  November. 

South  of  10°S  the  winds  a  directed  towards  WNW  with  strongest  winds 
in  September  and  October.  North  of  the  equator  the  monsoon  winds  blow  towards 
SW  from  November  to  March  with  maximum  wind  stress  in  January.  The  sum- 
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Figure  4 ■  Wind  stress  and  its  curl  over  the  Indian  Ocean  for  May.  The 

contour  interval  is  0.01  N/rr, t2  and  f.O-KT^  N/m^ . 


Figure  5.  Wind,  stress  and  its  curl  over  the  Indian  Ocean  for  August.  The 
contour  interval  is  0.01  N/rr?  and  4-0- 10~^  N/vr?. 
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mer  monsoon  is  much  more  intense,  with  northeastward  winds  from  May  through 
September.  A  bicubic  spline  was  used  to  interpolate  the  pseudo-stress  to  the  model 
grid  from  the  2°  x  2°  grid  of  the  original  data  set.  Linear  interpolation  was  applied 
to  compute  the  wind  stress  at  a  given  time  step  from  the  monthly  average  values. 
To  avoid  excessive  generation  of  internal  waves  when  the  model  is  started  from  rest, 
the  wind  stress  field  is  multiplied  by  the  function 

(1  -exp(-f/T0))  (2-11) 

where  tq  =  20  days,  during  the  first  year  of  the  spin  up.  This  wind  forcing  was 
chosen  to  make  comparison  easier  with  the  results  of  Woodberry  et  al.  (1989),  who 
used  the  same  wind  forcing  to  drive  their  model.  Further  details  are  given  in  their 
work. 

With  the  wind  forcing  used  here,  rapid  wind  changes  are  not  resolved. 
Using  scatterometer  data  with  6  days  resolution  compared  to  COADS  data  with 
30  days  resolution,  Perigaud  and  Delecluse  (1988)  found  that  the  response  of  a  1.5 
non-linear  reduced-gravity  model  was  insensitive  to  the  high  frequency  part  of  the 
forcing. 

The  model  is  spun  up  from  rest  by  applying  the  annual  wind  stress  cycle 
repeatedly,  until  a  quasi-periodic  solution  is  obtained.  From  experience  with  1.5 
layer  reduced-gravity  models  (Luther  and  O’Brien,  1985,  Luther  et  al.,  1985  and 
Woodberry,  et  al.  1989)  a  spin-up  time  of  less  than  20  years  can  be  expected  with 
three  baroclinic  modes.  In  the  two  first  of  these  studies  a  spin  up  time  of  three 
years  was  reported  to  be  adequate  for  the  Somali  Current.  This  was  due  to  the 
limited  area  of  the  model.  In  the  last  work,  the  spin  up  of  the  southern  part  of 
the  Indian  Ocean  required  approximately  7  years  of  integration  to  obtain  a  periodic 
solution.  Here,  since  higher  vertical  modes  propagate  slower,  the  spin  up  time  can 
be  expected  to  be  longer. 


3.  RESULTS 

Before  we  focus  on  the  Somali  Current,  the  ability  of  the  model  to  simulate 
the  general  circulation  of  the  Indian  Ocean  is  briefly  discussed.  The  results  described 
here  are  based  on  the  layer  thicknesses  and  velocity  fields  at  day  16  each  month, 
averaged  from  year  11  through  year  20  of  the  spin  up.  The  equatorial  currents  and 
counter  currents  as  well  as  the  upwelling  or  cooling  of  the  Arabian  Sea  west  of  India 
are  reproduced  well  by  the  model.  Climatological  maps  of  the  surface  currents  can 
be  found  in  VVyrtki  (1971),  Knox  and  Anderson  (1985),  and  Knox  (1987). 

3.1.  General  Circulation  of  the  Indian  Ocean 

The  souou  .quatorial  current  (SEC)  is  the  most  persistent  current  in  the 

model  as  well  ds  in  the  Indian  Ocean.  Its  westward  flow  reach  from  about  10°  S  to 
20°  S,  a  couple  of  degrees  further  south  during  the  northern  summer.  This  widening 
is  iu  agreement  with  observations,  which  also  show  an  increase  in  transport  from 
33  Sv  to  39  Sv  during  the  same  period  (e.g.,  Wyrtki,  1971,  Schott,  1983).  The 
average  transport  in  the  model  is  25.7  Sv,  while  the  1.5  layer  model  has  an  average 
transport  of  25.0  Sv.  While  the  single  layer  solution  has  a  small  seasonal  signal, 
the  transport  varies  from  19  Sv  in  April  to  31.6  Sv  in  September  for  the  3.5  layer 
case.  For  the  upper  layer  the  transports  for  the  same  months  are  21.0  Sv  and  26.4 
Sv,  respectively,  but  the  maximum  is  seen  in  October  and  November  with  26.7  Sv. 

The  SEC  separates  at  20°S  into  a  southward  and  northward  branch,  the 
latter  feeding  the  northward  flowing  East  African  Coastal  Current  (EAC),  which 
is  very  important  for  the  Somali  Current  by  either  enhancing  or  opposing  its  flow. 
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Figure  7.  Current  velocity  in  the  upper  layer  of  the  Indian  Ocean.  Ten  year 
average  of  February  16.  top:  1.5  layer  model,  bottom:  3.5  layer  model 


Figure  8.  Current  velocity  in  the  second  (top)  and  third  ( bottom )  layer  of  the 
Indian  Ocean.  Ten  year  average  of  February  16. 
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Figure  9.  Ten  year  average  February  16.  Depth  to  the  bottom  of  layer  1 
(a),  layer  2  (h),  layer  3 ,  (c)  and  surface  elevation  (d)  in  the  3.5  layer  model. 
Panels  (e)  and  (f)  show  the  surface  elevation  and  upper  layer  thickness  in  the 
1.5  layer  model. 
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Figure  10.  Current  velocity  in  the  upper  layer  of  the  Indian  Ocean.  Ten  year 
average  of  May  16.  top:  1.5  layer  model,  bottom:  3.5  layer  model 
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Figure  12.  Ten  year  average  May  16.  Depth  to  the  bottom  of  layer  1  (a),  layer 
2  (b)  layer  3,  (c)  and  surface  elevation  (d)  in  the  3.5  layer  model.  Panels  (e) 
and  ([)  show  the  surface  elevation  and  upper  layer  thickness  in  the  1.5  layer 
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Figure  13.  Current  velocity  in  the  upper  layer  of  the  Indian  Ocean.  Ten  year 
average  of  August  16.  top:  1.5  layer  model,  bottom:  3.5  layer  model 
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Figure  14-  Current  velocity  in  the  second  (top)  and  third  ( bottom )  layer  of 
the  Indian  Ocean.  Ten  year  average  of  August  16. 
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Figure,  15.  Ten  year  average  August  16.  Depth  to  the  bottom  of  layer  1 
(a),  layer  2  ( b ),  layer  3,  (c)  and  surface  elevation  (d)  in  the  3.5  layer  model. 
Panels  (e)  and  (f)  show  the  surface  elevation  and  upper  layer  thickness  in  the 
1.5  layer  model. 
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Figure  16.  Current  velocity  in  the  upper  layer  of  the  Indian  Ocean.  Ten  year 
average  of  November  16.  top:  1.5  layer  model,  bottom:  3.5  layer  model 
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Figure  17.  Current  velocity  in  the  second  (top)  and  third  ( bottom )  layer  of 
the  Indian  Ocean.  Ten  year  average  of  November  16. 
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e.g.,  Anderson  and  Moore,  (1979).  Its  average  transport  at  8°S  is  17  Sv  in  our 
models.  The  relative  strength  of  the  southward  branch  compared  to  the  northward 
along  the  east  coast  of  Madagascar  is  determined  in  part  by  a  parameter  in  the 
open  boundary  condition.  It  was  chosen  in  order  to  obtain  a  southward  flow  in 
agreement  with  observations.  The  average  volume  transport  is  12  Sv  in  the  models. 
A  second  branching  is  enforced  at  the  coast  of  East  Africa,  where  the  strength  of 
the  southward  flowing  Mozambique  Current  relative  to  the  East  African  Current 
also  depend  on  a  parameter  choice  in  the  open  boundary  condition.  The  southward 
transport  is  about  5  Sv.  The  open  boundaries  and  the  associated  parameters  is 
discussed  in  Appendix  B. 

During  the  winter  the  northern  hemishere  equivalent  of  the  SEC,  the  north 
equatorial  current  (NEC)  is  present  from  5°N  to  the  equator,  although  it  is  a  much 
weaker  current.  During  the  spring,  after  the  first  onset  of  the  southwest  monsoon 
in  April,  this  current  becomes  erratic,  and  joins  the  Equatorial  Counter  Current 
(EEC)  to  the  south  of  it,  to  form  the  eastward  flowing  Southeast  Monsoon  Current. 
This  reaches  from  3°N  to  9°S,  being  strongest  and  most  persistent  south  of  the 
equator,  in  the  region  2°S  -  7°S.  where  the  EEC  flows  during  the  winter.  The 
main  source  of  the  EEC,  however,  is  the  EAC  which  turns  offshore  just  south  of 
the  equator  during  the  northeast  monsoon,  and  north  of  the  equator  during  the 
northern  summer.  In  October  after  the  transition  to  the  northeast  monsoon,  the 
NEC  returns  to  a  westward  flow.  The  currents  in  each  layer,  the  depths  of  the 
isopycnals  and  the  surface  elevation  computed  by  (2-9)  for  the  two  models  are 
shown  as  in  Fig. 7- 18  for  February,  May,  August  and  November.  The  agreement  for 
the  upper  layer  model  results  and  the  dynamic  height  at  100  m  relative  to  1000  m 
in  Wyrtki  (1971)  is  good. 
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The  currents  in  the  second  and  third  layers  arc  in  the  same  direction 
along  the  equator,  but  generally  in  opposite  direction  to  the  upper  layer  flow,  e.g., 
second  vertical  mode.  These  deep  currents  change  direction  four  times  during  the 
year  (Figs.  8,11,14  and  17).  At  the  eastern  part  of  the  ocean  westward  pulses  of 
undercurrents  appear  in  April  and  November,  while  eastward  undercurrents  are 
generated  there  in  January  and  July.  The  patch  of  strong  undercurrents  move 
westward  with  a  phase  speed  of  about  40  cm/s,  which  is  in  agreement  with  linear 
theory  for  the  first  horizontal  Rossby  mode.  The  reversal  occur  first  in  the  lowest 
layer  indicating  an  upward  slope  to  the  east  of  lines  of  constant  phase,  e.g.,  upward 
phase  propagation. 

The  model  currents  reverse  with  depth  as  in  the  observations  of  Luyten 
and  Swallow  (1976)  who  found  an  eastward  upper  layer  flow  and  westward  flow 
below  in  May  and  June  1976.  The  semi-annual  changes  in  response  to  the  Monsoon 
wind  below  the  thermocline  at  750  m,  reported  by  Luyten  and  Roemmich  (1982), 
and  reproduced  here  as  Fig.  19,  are  in  excellent  agreement  with  the  model  response 
in  the  second  and  third  layer.  They  also  found  westward  and  upward  propagation 
of  phase.  Gent,  et  al.,  (1983)  presented  an  analytical  linear  model,  which  forced 
by  the  semi-annual  zonal  component  of  the  Hellerman  and  Rosenstein  wind  stress 
reproduced  these  reversals.  They  found  that  reflection  at  lateral  walls  was  essential 
for  a  realistic  result. 

Along  the  coastline  in  the  eastern  part  of  the  Indian  Ocean  we  have  high 
pressure  in  the  surface  during  the  summer  and  a  low  during  the  winter.  During 
the  summer  the  high  along  the  north  west  coast  of  Australia  drives  the  southward 
flowing  coastal  Leeuwin  Current.  West  of  India,  in  the  Bay  of  Bengal,  a  clockwise 
gyre  is  present  in  the  late  winter  and  early  spring,  while  an  anticlockwise  circulation 
is  found  in  late  summer  and  early  fall  in  agreement  with  observed  flow  patterns. 


40 


(5  OAY  AVERAGE  CURRENTS 
750  METERS 


.  -OS 

i 

-r-rj, - 

, ; ,  -  - - 

~?  r.  i  \ - 

«  ~ 

\  — 

;*  a 

-Tr-TTT - 

I  -31 

X  i:  -  * 

'  -  -  ' 

4 

rr  ,-:'i - 

!  ' 

' '  ~~  V  ^ 

l\  / 

\ 

,v 

1 

/  —  __r-  _ 

1 

-  £  -  -C  " 

^  ^ 

•  ! - 

.  -a 

.1  -14 

-  i  -  '>  • 

'  ^  <  • 

. . .  •  1 1 

.  ... 

i . .  - jj 

p 

IV  -,4 

0  WO  0  40 

««  Cn/Jte 


Figure  19.  Time  series  of  15  day  average  equatorial  current  at  750  m.  Data 
is  from  ten  moorings  from  f7°E  to  61° E  covering  the  time  period  May  1979 
to  April  1980.  From  Luyten  and  Roemmich  (1982). 
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In  the  Arabian  Sea,  west  of  India,  a  high  pressure  develops  during  the 
spring,  creating  an  anticyclonic  gyre  in  the  western  part  of  that  area.  In  May 
upwelling  and  an  associated  southward  current  starts  along  inc  west  coast  of  India. 
In  June  coastal  upwelling  is  also  found  along  the  Somali  Coast,  and  the  shallowing 
of  the  thermocline  is  seen  in  the  model  along  the  coasts  and  the  northern  part  of  the 
Arabian  Sea  from  July  to  October.  This  low  surrounds  a  a  clockwise  gyre,  which  is 
found  east  of  the  Island  of  Socotra.  It  is  fed  from  the  south  of  the  northward  flowing 
Somali  Current  with  an  offshore  return  flow.  These  model  results  (Fig. 9, 12, 15  and 
IS)  show  in  general  the  same  pattern  as  the  climatological  thermocline  depths  given 
by  Molinari,  Festa  and  Swallow,  (1986).  Their  maps  for  February  and  August  are 
shown  as  Fig.  20. 

The  variability  due  to  wind  forcing  is  clearly  seen  from  the  figures  showing 
the  total  depth  of  the  moving  layers.  Westward  propagating  planetary  waves  are 
radiated  from  the  eastern  boundaries  into  the  interior  of  the  ocean.  For  instance, 
annual  Rossby  waves  are  emitted  from  the  North  Coast  of  Australia,  the  Java 
coast  and  the  west  coast  of  India.  The  generation  of  highs  in  the  winter  and  lows 
during  the  summer  along  the  Indian  west  coast  is  the  opposite  of  what  would  be 
expected  from  local  coastal  upwelling.  The  changes  in  upper  layer  thickness  is  more 
likely  associated  with  the  changes  in  wind  stress  curl  seen  over  Sri  Lanka.  The 
pressure  pertubation  generated  there  will  propagate  northward  as  coastal  Kelvin 
waves  with  westward  radiation  of  planetary  waves,  (see  for  instance,  Gill,  1982).  At 
the  western  boundary  reflection  causes  the  wave  length  to  decrease  by  a  large  factor 
(e.g.,  Pedlosky,  1979),  which  explains  the  small  length  scale  and  eddy  activity  seen 
along  that  boundary.  Partial  reflection  takes  place  where  islands  are  included  in 
the  model. 
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At  the  equator  near  the  western  boundary  we  find  the  highest  variability 
in  the  model  solution.  The  solution  show  strong  oscillations  from  July  to  February 
in  the  meridional  velocity  component,  (Fig.  21). 

From  the  figure  we  find  a  period  of  2S  days  and  a  westward  phase  velocity 
of  37  cm/s.  corresponding  to  a  wave  length  of  900  km.  The  group  velocity  is  about  30 
cm/s  and  eastward.  The  phase  in  the  third  layer  leads  the  second  layer  which  implies 
upward  phase  propagation,  and  accordingly,  downward  energy  propagation.  Using 
these  values  in  the  linear  dispersion  relation  for  equatorial  waves,  we  can  identify 
the  oscillation  as  mixed  planetary-gravity  waves  (often  referred  to  as  Yanai  waves). 
Observations  show  that  such  oscillations  are  found  in  all  oceans,  (see  Weisberg. 
1987,  for  a  review).  They  were  first  observed  in  the  Indian  Ocean  by  Luyten  and 
Roemmich,  (1982).  Kindle  and  Thompson,  (1989)  studied  the  Yanai  waves  in  detail 
using  their  1.5  layer  reduced-gravity  model  of  the  Indian  Ocean.  They  forced  the 
model  by  the  original  Hellerman  and  Rosenstein  wind  stress  and  found  periods  in 
the  range  20-30  days  and  wave  lengths  from  800  km  to  1400  km.  The  generation 
was  attributed  to  barotropic  instability  of  the  Southern  Gyre  during  the  summer 
monsoor.  and  associated  with  meandering  of  the  East  African  Current  during  the 
winter.  Maps  of  potential  vorticity  for  year  21,  shown  in  section  3.4,  show  that  the 
nessessary  condition  for  this  type  of  instability  is  satisfied,  so  the  same  mechanism 
is  suggested  to  be  responsible  in  the  multiple  layer  case. 

The  model  solution  is  not  strictly  periodic.  A  comparison  of  individual 
years  shown  that  eddies  and  currents  appear  at  slightly  different  times  and  positions, 
although  the  overall  pattern  is  the  same.  This  internal  variability  is  illustrated  for 
February  16  and  August  16,  Fig.  22  and  Fig.  23,  by  the  variance  of  the  velocity 
and  standard  deviation  of  layer  thickness  anomalies  from  the  average  fields  at  these 
dates  taken  over  the  last  ten  years  of  integration.  The  first  quantity  can  be  regarded 
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Figure  22.  Based  on  ten  year  simulation  of  February  16:  log\Q  of  the  velocity 
components  (m/s)  in  layer  l  to  3,  ( panel  a-c).  Standard  deviation  of  layer 
thickness  1  to  3,  in  meters.  ( panel  d~f). 


Figure  23.  Based  on  ten  year  simulation  of  August  16:  log\{)  of  the  velocity 
components  (m/s)  in  layer  1  to  3,  (panel  a-c).  Standard  deviation  of  layer 
thickness  1  to  3,  in  meters,  (panel  d-f). 


as  proportional  to  an  eddy  kinetic  energy,  while  the  latter  is  related  to  the  square 
root  of  the  eddy  potential  energy.  The  variability  in  the  western  equatorial  region 
is  highest  in  the  winter  and  is  caused  by  the  Yanai  waves  discussed  above.  In  the 
Somali  current  we  find  that  the  solution  is  nearly  periodic.  The  variation  in  upper 
layer  thickness  show  the  same  pattern  found  by  Luther  and  O’Brien,  (1989)  for  their 
1.5  layer  model.  Along  the  east  coast  of  Madagascar  we  find  high  values  in  the  layer 
thickness  deviation,  in  particular  in  the  lower  layers.  The  natural  variability  may 
not  be  as  large  as  shown,  if  the  model  is  not  fully  spun  up  after  11  years  in  the 
southernmost  part  of  the  basin.  However,  the  recent  numerical  calculation  of  the 
world  ocean  circulation  by  Semtner  and  Chervin,  (1988),  using  the  annual  mean  of 
the  Hellerman-  Rosenstein  wind  stress,  exhibits  strong  eddy  activity  in  that  region. 

The  upper  layer  solution  of  the  1.5  layer  model  is  not  much  different  from 
the  3.5  layer  solution.  A  comparison  of  the  surface  elevations,  depth  of  the  upper 
layer  thickness  and  the  upper  layer  currents  show  little  differences  in  the  large  scale 
circulation.  Since  the  initial  value  of  the  phase  speed  of  the  first  baroclinic  mode 
differs  by  about  20%,  an  the  very  long  integration  taken  into  consideration,  this 
result  is  not  obvious.  Further  it  should  be  noted  that  the  results  from  this  1.5  layer 
model  are  nearly  identical  to  those  presented  by  Woodberry,  et  al.,  (1989)  although 
frictional  terms  and  open  boundary  conditions  are  implemented  differently. 

3.2.  Upper  layer  flow  of  the  Somali  Current 

Here  the  results  of  both  the  1.5  layer  model  and  the  3.5  layer  model  will 

be  presented.  An  average  of  the  flow  from  year  10  to  year  20  of  the  spin  up  will 
be  discussed  in  this  section.  The  discussion  of  the  result  are  based  on  the  3.5  layer 
version,  but  differences  in  the  solution  of  the  single  layer  model  are  pointed  out. 

The  model  produces  a  southward  Somali  Current  from  November  through 
March.  The  current  is  approximately  150  km  wide,  with  a  typical  alongshore  speed 


48 


of  20  cm/s,  and  extends  as  a  continuous  current  from  10°N  to  the  equator  in  that 
period.  During  January  to  March  the  current  is  southward  to  2°S,  before  it  meets 
the  northward  flowing  East  African  Current  and  turns  offshore  into  the  equatorial 
counter  current.  The  maximal  southward  transport,  14.4  Sv,  occurs  in  the  end  of 
January  in  the  3.5  layer  case.  With  one  layer  the  maximal  transport  is  8.6  Sv  in 
early  December. 

The  Somali  Current  is  fed  by  on  shore  flow  along  10°N  and  an  inflow- 
perpendicular  to  the  coast  at  5°N.  These  model  results  for  the  winter  Somali  Current 
match  the  observations  as  given  in  the  introduction,  or  as  presented  in  the  dynamical 
height  charts  by  Wyrtki  (1971).  The  single  layer  model  has  a  southward  flow 
only  to  the  beginning  of  March,  and  an  offshore  northward  flow  not  found  in  the 
climatological  charts.  The  reason  for  this  difference  in  models  behavior  is  due  to 
differences  in  response  to  the  decay  of  the  summer  monsoon  and  will  be  discussed 
later. 

During  the  northeast  monsoon,  the  model  has  a  decaying  cyclonic  eddy 
north  of  Socotra,  while  a  weak  anticyclonic  gyre  just  south  of  the  island  intensifies. 
With  one  layer  the  circulations  are  in  the  opposite  directions  (Fig.  24).  Judging 
from  thermocline  depth  observations  of  Molinari,  et  al.,  (1986a, 1986b),  the  multi 
layer  result  is  the  most  realistic.  The  annual  variation  of  the  thermocline  depth 
exhibits  a  maximum  in  March  and  a  minimum  in  September  with  -  amplitude 
of  35  m.  The  model  calculations  show  an  amplitude  of  30  m,  but  the  maximum 
thickness  occurs  in  January  in  the  1.5  layer  case. 

During  the  transition  period  to  the  southwest  monsoon,  from  March  to 
early  May,  the  multilayer  model  results  are  not  in  as  good  agreement  with  the 
observations.  The  weak  northward  surface  flow,  north  of  5°N  in  March,  and  along 
the  entire  coast  in  April  as  shown  in  climatological  maps  and  corroborated  by 
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Schott  and  Quadfasel  (1982)  with  current  measurements  during  1979,  is  missing 
in  that  version  of  the  model.  Schott  (1983)  attributes  this  to  “strong  wind  stress 
curl  distribution  off  Somalia,’'  but  the  wind  stress  used  to  force  the  model  show 
no  evidence  of  such  curl  and  one  can  expect  that  the  lack  of  this  early  onset  may 
be  due  to  insufficient  spatial  and  temporal  resolution  of  the  winds.  However,  the 
single  layer  version  does  have  northwards  currents  north  of  5°  at  that  time.  The 
early  onset  of  the  summer  monsoon  in  April,  due  to  local  northward  winds  south 
of  the  equator  as  reported  by  Leetmaa  (1972,1973)  is  seen  in  the  model  solutions. 

In  May  and  June  coastal  upwelling  starts  a  few  degrees  north  of  the  equator 
and  is  later  found  along  all  of  the  Arabian  Peninsula.  The  currents  are  to  the  north 
along  the  coast.  Since  the  local  wind  stress  is  strong  at  this  time,  we  can  expect 
the  models  to  give  similar  results  as  seen  in  Fig.  25.  The  Great  Whirl  forms  and 
migrates  to  the  north  during  July  and  August  and  the  anticyclonic  Socotra  eddy 
forms  east  of  that  island.  The  observations  show  that  the  Socotra  eddy  forms 
regularly  each  year  at  this  time,  and  that  its  high  salinity  distinguish  it  from  the 
Great  Whirl  (Bruce,  1979).  In  the  model  solution  the  eddy  forms  between  the 
northeastern  branch  of  the  Somali  current  and  a  return  flow  from  the  north.  Its 
potential  vorticity  corresponds  to  a  higher  latitude  than  that  of  the  Great  Whirl, 
which  suggest  that  it  consists  of  different  water  mass.  The  southern  gyre,  with 
offshore  flow  just  north  of  the  equator  (Fig.  26)  is  seen  in  the  solution  during 
the  entire  summer.  At  the  height  of  the  Monsoon  the  offshore  flow  is  at  2.5°N  in 
agreement  with  observations.  Fig.  27  shows  the  detailed  observations  of  the  surface 
currents  and  salinity  from  1979. 

In  most  years  a  northward  migration  of  the  southern  gyre  followed  by  a 
merging  with  the  Great  Whirl  is  observed  (Swallow  and  Fieux,  1982).  This  does 
not  happen  in  this  simulation.  Using  a  1.5  layer  model  very  similar  to  the  one  used 


Figure  24 .  Ten  year  average  February  16.  Velocity  and  depth  to  the  bottom 
of  layer  1,  (a.),  layer  2,  ( b ),  layer  3,  (c)  and  surface  elevation,  (d)  in  the  3.5 
layer  model.  Panels  ( e )  and  (f)  show  the  surface  elevation  and  upper  layer 
thickness  with  velocity  vectors  in  the  1.5  layer  model. 


Figure  25.  Ten  year  average  May  16.  Velocity  and  depth  to  the  bottom  of  layer 
1,  (a),  layer  2,  (b),  layer  3,  (c)  and  surface  elevation,  (d)  in  the  3.5  layer 
model.  Panels  (e)  and  (f)  show  the  surface  elevation  and  upper  layer  thickness 
with  velocity  vectors  in  the  1.5  layer  model. 
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Figure  26.  Ten  year  average  August  16.  Velocity  and  depth  to  the  bottom  of 
layer  1,  (a.),  layer  2,  (b),  layer  3,  (c)  and  surface  elevation,  ( d )  in  the  3.5 
layer  model.  Panels  (e)  and  (f)  show  the  surface  elevation  and  upper  layer 
thickness  with  velocity  vectors  in  the  1.5  layer  model. 
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Figure  21.  Observed  surface  current  vectors,  salinities  and  a  drifter  trajectory 
during  June- July  1979.  Surface  winds  are  shown  in  right  panel.  From  Diiing, 
Molinari  and  Swallow  (1980). 
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here.  Luther  and  O'Brien  (1989)  produced  this  merging  when  applied  to  observed 
monthly  wind  for  most  years.  However,  using  the  same  model  merging  did  not  occur 
when  the  Hellerman  and  Rosenstein  wind  stress  was  applied  as  forcing,  (Luther, 
personal  communication).  The  results  here,  in  spite  climatological  averaged  winds 
are  used,  does  therefore  correspond  to  the  less  common  situation  where  the  southern 
gyre  stays  near  the  equator.  During  the  summer  the  transport  in  the  model  reaches 
a  maximum  northeastward  transport  of  19.2  Sv  in  July  with  three  layers.  The  1.5 
layer  version  has  a  maximum  of  18.3  Sv  in  August. 

In  September  a  cyclonic  eddy  is  created  between  the  Great  Whirl  and  the 
Socotra  Eddy,  after  a  wedge  of  cold  water  is  advected  offshore  by  the  Great  Whirl. 
At  this  time  the  two  versions  of  the  model  starts  to  show  different  results.  In  the 
single  layer  case  a  small  anticyclonic  eddy  forms  north  west  of  Socotra.  In  October 
the  3.5  layer  version  still  has  a  northward  Somali  current,  and  the  Great  Whirl 
has  weakened  considerably.  With  one  active  layer,  the  whirl  is  stronger,  positioned 
further  north  and  northward  flow  is  confined  to  north  of  5°N. 

In  November  the  northeast  monsoon  has  arrived,  but  a  northward  coastal 
flow  is  still  found  north  of  7°N  as  indicated  by  observations.  In  the  multilayer 
version,  the  anticyclonic  eddies  have  essentially  disappeared,  while  they  are  still 
present  in  the  single  layer  model  (Fig. 28). 

The  main  reason  for  these  differences  between  the  1.5  and  3.5  layer  model 
during  the  transition  from  summer  to  winter  monsoon,  is  that  kinetic  energy  can 
be  transferred  vertically  in  the  latter  version.  During  this  time  we  see  that  the 
currents  in  the  second  and  third  layer  are  strongly  accelerated. 

3.3.  The  undercurrents  off  Somalia 

The  observations  of  the  subsurface  flow  are  sparse.  Adding  to  this  the  fact 
that  the  horizontal  correlation  scale  decrease  from  160  km  at  the  surface  to  90  km 
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at  700  m,  (Leetmaa,  et  al.,  1980),  our  knowledge  of  the  actual  undercurrent  are 
unfortunately  rather  limited.  Furthermore,  observations  from  shipboard  profiling 
obtained  in  May  1979,  differs  from  the  long  term  measurements  composed  of  current 
meter  records  near  5°N  (Quadfasel  and  Schott,  1983,  Fig.  29). 

Recent  results,  Schott  (1987),  from  the  French-U.S.  cooperative  program 
in  the  western  Indian  Ocean  show  that  the  undercurrents  at  the  equator  differs  from 
these  current  meter  results.  The  model  results  in  the  lower  layers  are  displayed  in 
Figs.  24,  25,  26,  and  28. 

The  narrow  undercurrents  require  very  fine  resolution  of  the  numerical 
model,  and  no  ocean  model  is  available  which  can  reproduce  all  observations.  This 
is  also  true  for  the  present  study,  but  as  will  be  clear  from  the  discussion  below, 
the  model  results  are  in  agreement  with  many  of  the  observed  characteristics  of  the 
undercurrents. 

During  the  winter  monsoon,  from  the  onset  in  November  until  it  weakens 
in  February,  a  subsurface  cyclonic  gyre  in  the  model  solution  is  causing  a  southward 
current  along  the  Somali  Coast.  It  reaches  from  from  10°N  to  2°N  in  November  and 
is  weakening  with  its  northern  latitude  decreasing  to  7°N  in  mid  February.  The  flow 
turns  offshore  and  returns  north  about  200  km  away  from  the  coast.  At  the  equator 
the  deeper  currents  are  northward  except  for  December.  In  January  and  February 
the  northward  flow  is  very  weak.  Quadfasel  and  Schott’s  (1983)  results  from  5°N 
show  a  continuous  southward  flow.  The  model  flow  does  not  quite  reach  that  far 
south  during  the  last  month.  On  the  other  hand  the  rather  strong  southward  flow 
during  December  in  the  model’s  deeper  layers  is  in  contrast  to  the  observations. 

A  remarkable  agreement  is  found  between  the  structure  of  the  model  cur¬ 
rents  for  February  and  the  measurements  by  Bruce  and  Volkmann  (1969).  They 
found  a  deep  cyclonic  subsurface  eddy  with  maximum  velocities  at  52.5°E  (to  the 
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Figure  29.  Observed  alongshore  currents  at  5°  N  based  on  all  available  complete 
monthly  values  from  1975  to  1979.  (a.):  Vertical  profiles,  (b):  Time  senes  at 
two  levels.  From  Quadfasel  and  Schott,  (1983). 
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Figure  30.  Deep  eddy  observed  in  February  to  March  1965  along  P°  N. 
Geostrophic  currents  are  based  on  hydrography  with  reference  level  obtained 
from  neutrality-buoyant  floats.  From  Bruce  and  Volkmann,  (1969). 
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north)  and  at  55°  (to  the  south),  with  a  weaker  southward  current  (Fig.  30).  This 
is  position  where  the  Great  Whirl  is  found  in  the  late  summer.  It  has  therefore 
been  suggested,  that  since  the  surface  currents  of  the  Great  Whirl  can  be  observed 
a  couple  of  months  after  the  summer  monsoon  weakens,  this  submerged  eddy  might 
be  the  remnants  of  the  Great  Whirl,  which  persisted  beneath  the  surface,  Knox 
(1987).  It  will  be  clear  from  the  model  results  shown  later,  that  such  a  scenario  is 
not  supported  by  the  present  calculations. 

As  discussed  in  the  previous  section,  the  upper  layer  currents  are  not 
simulated  very  well  in  March  and  April,  so  discrepancies  may  also  be  expected  in  the 
layers  below.  In  March  we  find  in  the  second  layer  a  northward  undercurrent  north 
of  5°X  and  south  of  1°N,  with  southward  flow  in  between,  caused  by  an  cyclonic 
eddy  connected  with  strong  cross  equatorial  flow.  An  equatorial  undercurrent  with 
eastward  velocities  of  10  cm/s  are  found  east  of  50°E  in  the  second  layer,  while 
southwestward  flow  is  found  at  7°N.  The  latter  turns  westward  in  April  and  the 
cyclonic  eddy  moves  a  couple  of  degrees  northward.  The  equatorial  undercurrent 
has  nearly  disappeared  at  this  time.  Schott  and  Quadfasel  (1980)  reported  mean 
current  reversals  along  the  coast  during  the  spring.  At  5°N  and  7°N  the  mean 
current  was  southwestward,  but  at  6°  a  northeastward  flow  was  observed.  The 
climatological  results  given  in  Quadfasel  and  Schott  (1983)  suggest  a  southwestward 
current  of  15  cm/s  at  5°N,  found  in  the  model  only  where  the  cyclonic  eddy  touches 
the  coast.  The  complicated  flowpattern  in  the  model  suggests  that  the  alongshore 
intermediate  flow  may  depend  very  much  on  the  location  along  the  coast. 

In  May  and  early  June  the  upper  layer  flow  is  northeastward,  and  an  un¬ 
dercurrent  is  found  south  of  6°,  but  not  north  of  that  latitude  where  an  onshore 
undercurrent  is  seen  under  offshore  surface  flow.  A  westward  equatorial  jet  develops 
in  May  and  last  through  July.  This  jet  was  observed  in  May  and  June  in  1976  and 
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during  INDEX  in  1979  (Leetmaa,  et  al.,  1980,  1982).  The  authors  of  the  second 
paper  observed  a  splitting  of  the  equatorial  undercurrent  into  a  northeastward  and 
southeastward  branch  before  it  reached  the  coast,  with  an  intensification  of  the  un¬ 
dercurrent  in  June.  The  model  show  the  same  flow  pattern.  The  disappearence  of 
the  model  coastal  undercurrent  in  June  is  possibly  caused  by  arrival  of  this  equa¬ 
torial  flow.  However,  the  increasing  northward  winds,  causing  coastal  upwelling, 
set  up  an  onshore  pressure  gradient,  which  help  decrease  the  southward  deep  flow. 
Quadfasel  and  Schott  (1983)  found  that  the  Somali  undercurrent  ceased  in  May. 
The  observations  in  1979  by  the  same  authors,  Schott  and  Quadfasel  (1980),  show 
a  reversal  from  southwestward  to  northeastward  currents  in  early  July. 

A  relatively  detailed  synoptic  map  of  the  flow  at  700  m  between  May  10 
and  June  4  1979  is  given  by  Leetmaa,  et  al.,  (1980)  and  reproduced  here  as  Fig.  31. 
This  should  be  compared  with  the  model  solution  in  Fig.  33  from  May  28  year  21 
of  the  integration.  We  notice  that  the  actual  currents  are  larger  than  in  the  model 
by  a  factor  of  3  to  4  or  more.  However,  the  directions  are  in  good  correspondance 
with  the  observed  flow.  There  is  offshore  flow  north  of  8°N,  onshore  south  of  that 
latitude  with  northward  coastal  flow  between.  We  have  onshore  flow  at  the  equator 
with  southward  currents  along  the  coast.  The  observed  flow  from  2°N  to  5°N  is 
not  well  resolved,  but  exhibits  a  spatial  structure  not  inconsistent  with  the  model 
results.  The  northward  undercurrent  in  the  model  is  also  in  agreement  with  the 
findings  by  Schott  (1987).  However,  the  branch  of  the  clockwise  gyre  that  produce 
this  flow  in  the  model  is  further  south  in  early  May. 

Since  deep  currents  are  relatively  weak,  linear  theory  should  be  adequate  to 
explain  the  results,  unless  of  course  the  deep  flow  is  caused  by  instabilities.  For  this 
reason  it  is  encouraging  that  the  spatial  pattern  and  direction  of  the  currents  in  the 
lower  layer  agree  with  the  observations.  To  pursue  this  further,  we  will  consider  the 


Figure  31.  Velocities  at  700  m.  Solid  arrows  represent  data  observed  between 
May  16  and  22,  1979.  Dashed  arrows  indicate  data  from  May  26  to  31,  1979. 
Temperature  was  measured  along  ship  track  (light  line).  From  Leetma a,  et  a  1., 
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Figure  32.  Linear  mode  solution.  Section  perpendicular  to  the  coast  at  5°  N 
(left  panel)  and  section  of  alongshore  flow  44  km  from  the  coast.  The  magnitude 
of  the  alongshore  wind  stress  is  shown  above  (right  panel).  Contour  interval  is 
10  cm/s.  Shaded  areas  indicate  southward  flow.  From  McCreary  and  Kundu, 
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linear  theory  of  McCreary  and  Kundu  (1985).  Using  separation  into  normal  modes 
they  solve  the  steady  problem  along  a  western  boundary  forced  by  local  alongshore 
winds  north  of  the  equator.  Their  solution  should  be  valid  during  the  initial  phase  of 
the  summer  monsoon  when  local  forcing  is  dominating  and  before  non-linear  effects 
become  important.  Fig. 32  show  their  solution.  A  southwestward  undercurrent  is 
found  only  to  the  south  of  the  wind  forcing.  They  found  that  onshore  deep  flow 
along  the  southern  edge  of  the  wind  field  separated  into  a  northward  and  southward 
branch  at  the  coast.  We  also  note  that  offshore  countercurrents  are  present  in  their 
solutions.  The  numerical  model  solution  from  May  28,  year  21  (see  Fig.  33,  next 
section)  show  the  same  features. 

During  July  and  August  the  Great  Whirl  moves  northward  and  deepens 
to  reach  all  three  layers  in  the  model.  There  is  no  undercurrent  north  of  the  equa¬ 
tor.  which  agrees  well  with  the  July  observations  of  Quadfasel  and  Schott  (1983). 
However,  a  southward  recirculation  develops  offshore  and  moves  westward,  which 
produces  a  southward  coastal  undercurrent  in  September.  The  observations  show  an 
onse'  in  August.  The  southward  undercurrent  reach  to  4°N  where  it  turns  offshore 
into  .in  eastward  equatorial  undercurrent. 

In  September  a  cyclonic  eddy  forms  between  the  Great  Whirl  and  the 
Socc  ra  eddy  and  the  flow  changes  rapidly.  It  will  be  clear  from  the  figures  shown 
in  tb  e  next  section,  that  the  eddies  are  tilted  with  depth,  the  warm  eddies  to  the 
nortl  west  and  the  cold  to  the  south.  In  October  the  Great  Whirl  has  collapsed  and 
the  deep  part  been  replaced  by  a  large  subsurface  cyclonic  eddy.  The  formation  of 
this  flow  will  be  further  discussed  in  the  next  section. 

3.4.  The  Great  Whirl 

The  results  after  20  years  of  spin  up  will  be  used  in  this  section  to  analyze 
the  generation  and  decay  of  the  Great  Whirl.  The  model  result  was  recorded  every 
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6  days  and  most  of  the  conclusions  are  based  upon  observing  a  computer  animation 
of  the  evolution,  and  may  not  be  obvious  from  the  figures  shown  here.  The  results 
from  the  1.5  layer  model  are  also  from  year  21  of  the  integration.  We  define  the 
following  diagnostic  quantities  to  be  used  in  the  subsequent  analysis: 

Vertical  component  of  the  relative  vorticity  in  layer  j: 

Cj  =  V  x  Vj  (3-1) 

where  vj  is  the  velocity  vector  in  layer  j  and  V  x  the  curl  operator  in  spherical 
coordinates. 

Vertical  component  of  the  absolute  vorticity  in  layer  j: 

Caj  =Cj+f  (3  -  2) 

where  /  is  the  Coriolis  parameter. 

Vertical  component  of  the  potential  vorticity  in  layer  j: 

IT,-  =  Caj'/ffj  (3  -  3) 

Kinetic  energy  per  area  unit  for  each  layer: 

Ekin,j  =  2 Pfij  '  Vj Hj  (3  —  4) 

which  is  the  vertical  integral  of  the  kinetic  energy  density. 

Available  potential  energy: 

1  A-l 

A  =  Tjdiipl  -  Pa)v 2  +  E  “  Pj)(HJ  -  H0j)2)  (3  -  5) 

J=1 

where  pa  is  the  density  of  air  and  rj  the  surface  elevation  and  Hqj  the  layer  thickness 
of  the  ocean  at  rest.  The  first  term  which  involve  the  free  surface  is  negligble 
compared  to  the  contribution  from  the  isopycnal  displacements  of  the  deeper  layers. 
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Different  levels  of  zero  available  potential  can  be  chosen.  Here  a  constant  have  been 
preferred  compared  to  a  background  field,  i.e.,  a  climatology,  which  may  vary  in 
time  and  space. 

The  evolution  of  the  flow  during  the  summer  monsoon  is  shown  in  Fig.  33- 
48.  A  velocity  vector  is  shown  for  each  layer  for  every  16  grid  points  in  the  model, 
i.e.,  0.8°  between  the  vectors.  These  plots  show  also  the  depth  to  the  bottom  of 
each  layer.  For  the  first  layer  this  is  simply  the  layer  thickness,  which  corresponds 
to  the  thermocline  depth  in  this  model.  For  the  second  layer  this  depth  is  the 
sum  of  the  thicknesses  of  layer  1  and  2,  while  the  3.  layer  depth  is  defined  as  the 
total  thickness  of  the  active  layers.  The  surface  elevation,  which  is  proportional 
to  the  dynamic  pressure  for  the  first  layer,  is  calculated  according  to  (2-9)  and 
shown  with  the  velocities  for  layer  1.  The  quantities  defined  above  are  shown  for 
the  same  dates  and  help  identify  dynamic  active  regions.  For  instance,  westward 
propagating  Rossby  waves  are  seen  in  the  relative  vorticity  maps  than  in  the  velocity 
and  pressure  fields. 

In  late  May  and  early  June,  Fig.  33  and  Fig.  34,  the  Somali  Current  is 
well  developed.  The  flow  split  in  a  northward  and  offshore  branch  at  6°N.  The 
latter  turns  southeastward,  but  no  closed  eddy  is  formed.  This  flow  is  associated 
with  a  local  surface  depression  to  the  north  of  that  latitude.  A  westward  equatorial 
undercurrent  dominates  the  flow  in  the  two  deep  layers.  The  available  potential 
energy  is  low  except  in  the  gyre  south  of  the  equator. 

By  June  16  it  is  clear  that  a  weak  eddy  has  formed  offshore  at  5°  N, 
(Fig.  35).  Using  a  three-dimensional  numerical  model  with  a  meridional  coast,  Cox 
(1976)  showed  that  mean  kinetic  to  eddy  kinetic  energy  transfer  took  place  during 
the  formation  of  the  Great  Whirl.  The  mechanism  was  identified  as  horizontal  shear 
instability.  Water  with  low  potential  vorticity  is  advected  to  the  north,  Fig.  36.  The 
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fluid  also  conserves  its  absolute  vorticity,  sice  vortex  stretching  is  weak.  Because 
of  the  non-zonal  character  of  the  Somali  Current  and  the  presence  of  the  coast, 
classical  barotropic  instability  theory,  i.e.  Kuo’s  (1949)  nessessary  condition  for 
instability,  is  not  valid.  However,  flows  with  an  angle  to  the  east-west  direction 
generally  are  less  stable  than  zonal  flows,  Gill,  et  al.,  (1974),  Pedlosky  (1979).  The 
more  general  nessessary  condition  for  mixed  barotropic-baroclinic  instability  is  that 
the  potential  vorticity  gradient  of  the  mean  flow  changes  sign,  (Pedlosky,  1979),  i.e. 
that 

v3n  =  0  (3  -  6) 

somewhere  in  the  mean  flow.  The  subscript  indicates  that  the  gradient  operator 
is  three-dimensional.  A  sign  change  in  Vn  along  isopycnals  indicates  a  possibility 
of  barotropic  instability,  Charney  and  Stern,  (1962),  Boudra,  et  al.,  (19SS).  This 
condition  is  fullfilled  prior  to  the  generation  of  the  whirl,  Fig.  34,  and  a  vertical 
reversal  of  Vn  nessessary  for  baroclinic  instability  is  also  found.  An  increase  in  n 
is  found  in  the  second  layer.  In  case  of  the  1.5  layer  model  the  generation  of  the 
Great  Whirl  looks  very  similar.  The  main  difference  is  that  the  coastal  upwelling 
is  less  intense  due  to  the  larger  density  difference  between  the  upper  and  second 
layer  than  in  the  multilayer  case.  This  causes  weaker  alongshore  currents  in  that 
model.  A  high  pressure  situated  southeast  of  Socotra  compensates  for  this  north 
of  10°N  and  along  Socotra  and  produces  a  stronger  alongshore  flow  there,  but  a 
weaker  offshore  flow  is  found  at  6°N.  This  may  be  the  reason  that  the  generation 
of  the  eddy  is  delayed  by  about  a  week  compared  to  the  3.5  layer  case.  The  fact 
that  the  eddy  formation  is  essentially  the  same  in  the  two  models  and  no  decrease 
in  available  potential  energy  was  observed  indicates  that  baroclinic  instability  does 
not  play  a  role.  At  the  time  of  generation  the  available  potential  energy,  A,  was 
increasing  along  the  coast  as  well  as  the  kinetic  energy  of  the  Somali  Current,  i.e.. 


Figure  3\.  May  28,  year  21,  3.5  layer  model.  Kinetic  energy  (kJ/rrr ,  shading: 
E fc,n  >40)  in  layer  1,  (a),  and  layer  (2),  ( b ),  (100  J/m? ,  shading:  >50). 

Relative  vorticity  in  layer  1  (10~^s~^,  shading:  |£|  >10),(c ),  and  layer  2 
(10~^s~^ ,  shading:  |£|  >30), (d).  Potential  vorticity  in  layer  1,  (e),  and  layer 
2,  (f),  (10~8m~ls~l ,  2  <  fl  <  4  is  shaded). 


Figure  35.  June  16,  year  21,  3.5  layer  model.  Velocity  arid  depth  (m)  to  the 
bottom  of  layer  1,  (a),  layer  2,  (b),  layer  3,  (c),  velocity  of  upper  layer  and 
surface  elevation  (cm),  (d).  Panel  (e):  Available  potential  energy,  (kj/m2). 


Figure  36.  June  16,  year  21,  3.5  layer  model.  Kinetic  energy  (kj/rr?,  shading: 
Ekin  >40)  in  layer  1,  (a),  and  layer  (2),  (b),  (100  J/m 2,  shading:  Ekin  >50). 
Relative  vorticity  in  layer  1  ( 10~6s~l ,  shading:  |(|  >10),(c),  and  layer  2 
(10~'s  shading:  |£|  >30),  (d).  Potential  vorticity  in  layer  1,  (e),  and  layer 
2,  (f),  ( ltrHm-ls-[ ,  2  <  n  <  4  is  shaded). 


I 


an  ./?.  July  I(>.  iji  nr  Jl.  J.r>  luytr  wodtl.  If  lonlij  mid  d(  pth  (m)  to  the 
o in  of  huii  r  I.  Oil.  Inyt  r  J.  (b).  Itnji  r  ■  /,  (c),  nlnrity  of  u]>p<  r  layer  and 
fan  I  If  ration  'em),  hi).  I’nntl  if):  Arrulablt  pottntial  nuryy,  ( Ic.I/tn "I. 


Figure  38.  July  16,  year  21,  3.5  layer  model.  Kinetic  energy  (kJ/rrr ,  shading: 
Kkin  >40)  in  layer  1.  (a),  and  layer  (2),  (b),  (100  J/m 2,  shading:  Ekux  >50). 
[{flat ire  vorticity  in  layer  1  (10~C\ ,s_1,  shading:  |C|  >10),(c),  and  layer  2 
[10~'s~^,  shading:  |C|  >30).  (d).  Potential  vorticity  in  layer  1.  (e).  and  layer 
2.  (f).  f  I  2  <  II  <  l  is  shaded). 
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the  current  was  still  in  a  state  of  acceleration.  The  local  wind  is  northeastward  and 
intensifying,  with  negative  curl  offshore,  but  has  too  large  scales  to  be  able  to  force 
the  whirl  directly.  The  eddy  forms  at  the  place  along  the  coast  where  the  kinetic 
energy  in  the  current  is  maximum  and  the  largest  gradient  in  relative  vorticity. 

The  Great  Whirl  intensifies  and  moves  northward  along  the  coast.  It  is 
seen  in  the  second  layer  motion  12  days  after  the  generation  and  after  IS  days  in 
the  third  layer.  Since  the  field  is  only  saved  from  the  computation  every  six  days,  a 
more  accurate  estimate  cannot  be  given.  During  July  the  northward  migration  and 
strengthening  of  the  whirl  continues,  Figs.  37-38.  A  band  of  positive  vorticity  along 
the  coast  show  the  width  of  the  frictional  boundary  layer.  The  discontinuation  of 
the  northward  flow  at  2°N  due  to  offshore  flow  at  that  latitude  is  clearly  seen  at  this 
time.  Positive  vorticity  is  advected  offshore  by  the  Great  Whirl  on  its  northern  side 
and  wraps  around  the  “Golf  club”  like  negative  vorticity  region.  The  numerical 
solution  by  Cox  (1979)  showed  the  same  structure  of  the  relative  vorticity  field. 
The  kinetic  energy  fields  at  the  same  date,  show  that  the  strongest  currents  are 
found  where  the  Somali  Current  turns  offshore  into  the  whirl,  while  currents  are 
much  weaker  on  the  southern  side  in  agreement  with  observations  (see  for  instance 
Swallow,  et  al.,  1983). 

Fig.  38  (c)  shows  that  the  potential  vorticity  for  the  upper  layer  in  the 
center  of  the  Great  Whirl  correspond  to  that  at  2-3°N.  Some  flux  of  positive  relative 
vorticity  can  be  expected  from  the  boundary  layer  and  bv  wind  stress  curl,  but  from 
the  time  when  the  whirl  is  formed  in  mid  June  to  late  August,  the  potential  vorticity 
of  its  center  is  changed  very  little.  This  implies  that  the  whirl  is  capable  of  carrying 
mass  during  its  northward  migration. 


The  Somali  Current  is  stable  until  the  beginning  of  August.  A  low  of  the 
westward  propagating  Rossby  waves  reach  the  eastern  edge  of  *  he  Great  Whirl  and 
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cyclonic  motion  is  amplified  as  seen  in  Figs.  39-42,  showing  the  field  variables  and 
derived  quantities  at  August  10  and  August  22.  The  cyclonic  spinup  starts  in  the 
lower  layers  first.  The  result  is  that  a  small  cold  eddy  pinches  off  the  wedge  north 
of  the  Great  Whirl,  intensifying  the  flow  along  the  southern  edge  of  the  whirl,  in 
particular  in  the  deep  layers.  This  is  the  process  that  restores  the  undercurrent 
under  the  Great  Whirl.  The  available  potential  energy  of  the  whirl  is  essentially 
unchanged,  while  we  observe  a  decrease  in  kinetic  energy  in  the  upper  layer  and  an 
increase  in  the  second  layer. 

The  available  potential  energy  of  the  Great  Whirl  reach  its  maximum  at 
September  4,  Fig.  43.  The  axis  of  the  whirl  i11  now  tilting  towards  the  northwest 
with  depth.  A  new  strong  cyclonic  eddy  is  generated  between  the  whirl  and  Socotra. 
To  the  south  of  the  whirl  and  to  the  north  of  the  cyclonic  eddy,  weaker  eddies 
exist.  The  two  large  eddies  starts  to  migrate,  initially  southeastward.  Then  the 
cyclonic  eddy  moves  southward  to  be  nearly  at  the  same  latitude  as  the  Great 
Whirl,  fie.,  a  clockwise  rotation  of  the  axis  through  the  eddy  centers,  and  the  pair 
moves  westward.  The  initial  movement  is  most  likely  due  to  mutual  advection.  A 
theoretical  study  of  eddies  by  Matsuura  and  Yamagata  (1982),  showed  a  similar 
initial  clockwise  rotation  of  an  eddy  dipole  with  the  cyclonic  eddy  to  the  north. 

During  the  process  the  Great  Whirl  is  losing  potential  energy,  while  kinetic 
energy  is  increased  in  the  lower  layers  east  of  the  whirl,  Fig.  43-46.  The  low  to  the 
southeast  of  the  whirl,  due  to  the  seasonal  Rossby  wave  field,  is  amplified  at  the 
same  time.  The  nessessary  condition  for  baroclinic  instability  that  the  vertical 
gradient  of  potential  vorticity  changes  sign,  (Pedlosky,  1979)  is  satisfied  in  the 
area.  This  evolution  suggest  that  the  eddy  pair  is  subject  to  a  barotropic-baroclinic 
instability.  Naturally,  other  mechanisms  may  play  a  role  in  the  complicated  flow 
in  the  area.  The  northeastward  wind  stress  starts  to  relax  in  the  model  after  mid 
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August  although  the  pattern  of  wind  stress  and  its  curl  does  not  start  to  change  until 
after  September  16.  Model  results,  discussed  below,  show  that  this  relaxation  is  not 
the  major  cause  for  change  of  the  flow.  There  is  an  equatorial  influence  on  the  flow 
in  the  region:  The  deep  flow  south  of  the  eddies  is  connected  with  the  deep  eastward 
equatorial  undercurrent;  the  potential  vorticity  map  from  September  4  also  show 
an  influx  of  water  originated  south  of  the  equator,  which  form  a  small  anticvclonic 
eddy  at  6°N.  This  development  in  the  flow,  with  the  collapse  of  the  Great  Whirl, 
is  seen  by  comparing  the  results  for  September  4  and  September  22.  As  previously 
noted,  the  Hellerman-Rosenstein  wind  stress  results  in  the  anomalous  situation 
that  the  Southern  Gyre  remains  in  the  vicinity  of  the  equator.  The  formation  of 
the  small  eddy  which  joins  the  Great  Whirl  and  adds  negative  vorticity  is  the  only 
indication  of  this  tendency.  In  the  two  deep  layers  westward  propagation  of  the 
eddy  field  causes  a  increased  northward  current  through  the  strait  between  Somalia 
and  Socotra  and  a  negative  vorticity  flux  from  the  whirl  along  the  coast  of  Socotra 
to  its  northern  shore. 

Fig.  47-48  from  October  4  shows  that  the  subsurface  part  of  the  Great 
Whirl  has  been  detached  and  a  southwestward  undercurrent  now  exist  beneath  it. 
The  flow  south  of  Socotra  consist  now  of  a  series  of  eddies  along  a  nearly  zonally 
oriented  axis,  which  moves  westward  towards  the  Somali  coast.  The  whirl  decays 
further  during  October,  while  the  Socotra  Warm  eddy  intensifies  until  the  middle 
of  the  month,  when  it  becomes  subject  to  an  instability  process  similar  to  that 
described  above  for  the  Great  Whirl.  The  Socotra  eddy  moves  south  of  the  island 
in  its  westward  movement,  and  by  the  end  of  November  the  anticvclonic  eddies  have 
dissappcared. 

The  effect  of  the  weakening  monsoon  winds  iri  August  and  September 
on  this  development  is  investigated  by  recomputing  the  solution  for  year  21,  but 
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extending  the  length  of  the  summer  monsoon.  In  the  first  experiment,  case  1, 
observed  winds  was  used  until  August  16  and  then  kept  constant.  The  wind  stress 
pattern  and  it  curl  is  nearly  identical  in  July  and  August,  but  the  magnitudes  are 
larger  in  July.  In  this  case  we  have  a  monsoon  which  reaches  a  maximum  and  then 
relaxes  to  be  a  weaker  forcing.  The  second,  case  2,  kept  the  winds  constant  from 
July  16,  so  that  Somali  Current  spins  up  under  maximal  forcing. 

Both  cases  show  an  evolution  of  the  flow  very  similar  to  the  standard 
solution,  case  0,  forced  by  seasonal  winds,  until  late  September.  The  generation  of 
deep  cyclonic  eddies  to  the  southeast  of  the  whirl  takes  place  in  the  same  fashion, 
with  a  decrease  in  available  potential  energy  of  the  Great  Whirl,  but  the  intensity  of 
the  eddies  is  smaller  than  before.  The  explanation  for  this  may  be  as  follows:  With 
the  weakening  of  the  monsoon,  case  0  and  case  1,  the  magnitude  of  the  eastward 
pressure  gradient  and  the  associated  deep  northward  flow  along  the  coast  of  Somalia 
decreases.  This  allows  a  westward  movement  of  the  deep  part  of  the  whirl,  causing 
the  tilt  towards  the  northwest  with  depth  of  the  Great  Whirl,  observed  prior  to  this 
event,  to  increase.  Accordingly,  the  vertical  shear  under  the  whirl  increases  with 
decreasing  southwest  winds.  Thus  vertical  shear  instability  is  enhanced  compared 
to  the  case  where  the  wind  continues  to  blow  at  full  strength.  In  case  0  the  cyclonic 
eddy  generated  in  August  moves  westward  and  reach  the  Somali  coast,  which  creates 
an  undercurrent.  In  case  1  and  2  the  deep  flow  is  only  slowed  down,  but  continues 
to  be  northward.  However,  the  second  larger  deep  cyclonic  eddy  creates  a  coastal 
undercurrent  in  late  October  (case  1). 

The  furth  development  in  the  cases  with  extended  summer  monsoon 
season  gives  an  interesting  result.  Instead  of  decaying  as  in  case  0,  the  Great  Whirl 
recover  its  loss  of  energy  during  a  clockwise  rotation.  By  the  end  of  October  the 
solution  looks  like  the  standard  case  in  early  September,  and  a  new  cyclonic  eddy 
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is  generated.  The  instability  is  thus  an  important  mechanism  by  which  the  Great 
Whirl  is  able  to  dissipate  its  excess  energy. 

The  evolution  in  the  1.5  layer  standard  case  with  seasonal  forcing  is  dif¬ 
ferent  than  described  above.  When  the  low  pressure  arrives  in  the  beginning  of 
August,  a  southwestward  current  appears  between  this  low  and  the  Great  Whirl, 
which  is  squeezed  slightly  northward  as  seen  in  Fig.  49.  showing  the  results  from 
August  10.  No  cyclonic  eddy  is  generated  to  the  south  in  this  case,  but  the  cyclonic 
flow  in  the  wedge  between  the  whirl  and  the  Socotra  eddy  is  weaker  in  this  case. 
The  Great  Whirl  and  the  cyclonic  to  the  northeast  are  increasing  in  strength  until 
early  September,  then  they  both  start  to  lose  potential  and  kinetic  energy.  The 
available  potential  energy  decrease  fastest  for  the  Great  Whirl,  but  at  a  slower  rate 
than  in  the  3.5  layer  case.  The  solution  for  September  16  is  representative  during 
this  time,  (Fig.  51).  The  Socotra  eddy  is  intensified,  indicating  eastward  energy 
propagation.  The  eddies  south  and  east  of  Socotra,  most  clearly  seen  in  the  relative 
vorticity  plot,  move  slowly  westward  after  this  time.  Negative  vorticity  from  the 
Great  Whirl  leaks  northward  as  in  the  3.5  layer  case  (Fig.  46),  but  the  Socotra  eddy 
moves  to  the  north  of  the  Island  in  this  case  instead  of  south  of  it  as  before. 
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Figure  39.  August  10,  year  21,  3.5  layer  model.  Velocity  and  depth  (in)  to  the 
bottom  of  layer  1,  (a),  layer  2,  (b),  layer  3,  ( c ),  velocity  of  upper  layer  and 
jurface  elevation  (cm),  (d).  Panel  (e):  Available  potential  energy ,  ( kJ/m "*). 
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Figure  40.  August  10,  year  21,  3.5  layer  model.  Kinetic  energy  (kJ/m* ,  shad¬ 
ing:  £fcin  >40)  in  layer  1,  (a),  and  layer  (2),  (b),  (100  J /nr? ,  shading: 
Eif. in  >50).  Relative  vorticity  in  layer  1  (10~^s~^,  shading:  |£|  >10),(c), 
and  layer  2  (10~"^s~\  shading:  |£|  >30), (d).  Potential  vorticity  in  layer  1, 
(e),  and  layer  2,  (f),  ,  2  <  fl  <  4  is  shaded). 


Figure  J,  1 .  August  22.  year  21,  5  r>  layer  model.  Velocity  and  depth  (m)  to  the 
bottom  of  layer  1,  (a),  layer  2,  (b).  layer  3,  (c),  velocity  of  upper  layer  and 
surface  elevation  (cm),  (d).  Panel  (e):  Available  potential  energy.  ( kJ/m *). 
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Figure  42.  August  22,  year  21,  3.5  i 
ing:  Ekin  >40)  in  layer  1,  (a), 
Eiiin  >50).  Relative  vorticity  in  l 
and  layer  2  (10~7s~\  shading:  |(| 
(e), and  layer  2,  (f),  (10~^m~^s~^, 
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layer  model.  Kinetic  energy  ( kJ/m ^ ,  shad- 
and  layer  (2),  ( b ),  (100  J/m*,  shading: 
layer  1  (I0~^s~^,  shading:  |C|  >10),(c), 
>30),  (d).  Potential  vorticity  in  layer  1, 
2  <  FI  <  4  is  shaded). 


Figure  44 ■  September  4,  year  21,  3.5  layer  model.  Kinetic  energy  ( kJ/m 2, 
shading:  Ekin  >40)  in  layer  1,  (a),  and  layer  (2),  (b),  (100  J/m 2,  shading: 
^ kin  >50).  Relative  vorticity  in  layer  1  (I0~&s~x,  shading:  |(|  >10), (c),  and 
layer  2  (10~^s  shading:  |C|  >30), (d).  Potential  vorticity  in  layer  1.  (e), 
and  layer  2,  (f),  (10~8  m~ 1 s~ 1 ,  2  <  U  <  4  is  shaded). 
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Figure  46.  September  22,  year  21,  3.5  layer  model.  Kinetic  energy  (kj/rr? , 
shading:  E £,n  >40)  in  layer  1,  (a),  and  layer  (2),  (b),  (100  J/m 2,  shading: 
Ekin  >50).  Relative  vorticity  in  layer  1  (I0~^s  shading:  |£|  >10), (c),  and 
layer  2  (10~7s~^,  shading:  |£|  >30), (d).  Potential  vorticity  in  layer  1,  (e), 
and  layer  2,  (f),  (10~^m~^ s~^ , 2  <  FI  <  4  is  shaded). 


Figure  4  7,  October  4,  year  21,  3.5  layer  model.  Velocity  and  depth  (m)  to  the 
bottom  of  layer  l,  (a),  layer  2,  (b),  layer  3,  (c),  velocity  of  upper  layer  and 
surface  elevation  (cm),  (d).  Panel  (e):  Available  potential  energy,  (kJ/rr?). 
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Figure  J8.  October  4,  year  21,  3.5  layer  model.  Kinetic  energy  (kJ/m? ,  shad¬ 
ing:  >40)  in  layer  1,  (a),  and  layer  (2),  ( b ),  (100  J/rr i2,  shading: 


Ekin  >50).  Relative  vorticity  in  layer  1  (10~^s  shading:  |£|  >10),  (c), 
and  layer  2  (!0~^s~'L,  shading:  |£|  >30), (d).  Potential  vorticity  in  layer  1, 
(e),  and  layer  2,  (f),  (10~^m~ls~^ ,  2  <  II  <  4  is  shaded). 
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Figure  49.  August  10,  year  21,  1.5  layer  model.  Velocity  and  depth  (m)  to 
the  bottom  of  active  layer,  (a)  and  surface  elevation,  (b).  Kinetic  energy 
( kJ/m ? ,  shading:  E >40),  (c).  Available  potential  energy,  (kJ/n?),  (d). 
Relative  vorticity  (10~6s~l,  shading:  |(|  >10),(e).  Potential  vorticity.  (f). 
(l(r8m-ls-l,2  <n  <4J  is  shaded. 


Figure  50.  September  22,  year  21,  1.5  layer  model.  Velocity  and  depth  (m) 
to  the  bottom  of  active  layer,  (a)  and  surface  elevatior,  (b).  Kinetic  energy 
(kJ/m? ,  shading:  >40),  (c).  Available  potential  energy,  (kj/m?),  (d). 

Relative  vorticity  (10~^s~l,  shading:  J£|  >10),(e).  Potential  vorticity,  (f), 
( I0~sm~ls~l ,  2  <  II  <  4  is  shaded). 


4.  SUMMARY 

In  the  introduction  the  present  observational  knowledge  and  various  mod¬ 
elling  efforts  in  the  past  was  reviewed.  Major  oceanographic  field  programs  in  the 
sixties  and  late  seventies  have  increased  our  knowledge  about  the  equatorial  and 
northwestern  Indian  Ocean,  in  particular  the  Somali  Current  system.  While  the 
near-surface  flow  is  well  documented  during  the  summer  monsoon,  observations  are 
sparse  during  the  northern  winter  and  below  the  thermocline.  However,  the  avail¬ 
able  data  show  that  the  deeper  flow  is  quite  energetic  and  consist  of  eddies  with 
small  horizontal  scale  near  the  Somali  Coast  (e.g.,  Leetmaa,  et  ah,  1980,  Quadfasel 
and  Schott,  1982.). 

Previous  models  have  fallen  into  one  of  two  categories.  The  first  type  is 
that  of  process  models  with  simple  geometries  and  forcing,  which  explain  the  physics 
involved  in  certain  problems.  An  example  is  the  generation  of  the  Somali  Current 
system  by  local  wind,  (Hurlburt  and  Thompson,  1976,  Lin  and  Hurlburt,  1981). 
The  second  type  is  models,  which  simulate  the  circulation  with  realistic  geometries 
and  forcing.  Examples  are  the  models  by  Cox  (1970),  Luther  and  O’Brien,  (1985, 
1989)  and  Woodberry,  et  ah,  (1989).  The  first  of  these  studies  had  inadequate 
horizontal  resolution  to  resolve  eddies,  while  the  latter  authors  had  only  one  active 
layer.  The  present  study  falls  in  the  second  category  and  is  the  the  first  Indian 
Ocean  isopycnal  model  to  be  eddy  resolving  with  vertical  structure. 

The  objectives  of  this  study  are:  1.  Verify  that  the  model  reproduces  the 
observed  features  of  the  Indian  Ocean  circulation  and  in  the  Somali  Current  system. 
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2.  Use  the  model  results  to  explain  the  observed  undercurrent  structure.  3.  Find  a 
possible  explanation  for  the  decay  of  the  Great  Whirl. 


Section  2  presents  the  new  general  multi-layer  model.  It  is  formulated  in 
spherical  coordinates  with  few  approximations.  To  prevent  ventilation  of  the  upper 
layer,  its  initial  thickness  has  been  chosen  as  200  m,  and  an  entrainment  scheme 
was  added,  (McCreary  and  Kundu,  1988).  A  reduced-gravity  formulation  in  which 
the  pressure  gradient  vanish  in  the  lowest  layer  filters  out  the  external  modes.  The 
computer  implementation  allows  any  number  of  active  layers  to  be  specified  within 
memory  limitations. 

The  model  geometry  is  realistic  with  open  boundaries  in  the  east  and  the 
south.  The  forcing  is  the  climatological  seasonal  wind  stress.  The  initial  condition 
and  the  open  boundary  conditions  are  determined  using  normal  vertical  modes. 
The  formulation  of  the  open  boundary  conditions  for  multi-layer  models,  based  on 
Camerlengo  and  O’Brien,  (1980),  is  new.  The  fine  horizontal  resolution  ensures 
that  the  internal  radius  of  deformation  for  the  three  baroclinic  modes  in  the  model 
are  resolved. 

In  section  3  the  results  from  numerical  simulations  with  one  and  three 
active  layers  are  presented.  The  general  circulation  in  the  upper  layer  is  essentially 
the  same  as  found  by  1.5  layer  models  (Woodberry,  et  ah,  1989)  and  is  in  good 
agreement  with  observations.  The  semi-annual  reversals  of  the  undercurrents  ob¬ 
served  by  Luyten  and  Roemmich,  (1982),  at  intermediate  depths  are  reproduced. 
In  agreement  with  their  observations,  the  phase  propagation  is  westward  and  up¬ 
ward.  Generation  of  planetary  waves  along  eastern  land  boundaries  are  important 
for  the  interior  time-dependent  flow  and  the  solution  near  the  western  boundary. 
Here  the  long  planetary  waves  reflect  by  decreasing  their  wave  length  and  eddies 
and  meanders  are  generated  through  nonlinearities.  Compared  to  1.5  layer  models, 
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e.g.,  Luther  and  O’Brien,  (1989),  we  find  an  increased  variability  when  the  lower 
layers  are  active. 

At  the  equator  and  along  the  Somali  Coast  the  most  energetic  eddies 
are  found.  Yanai  waves  are  responsible  for  eastward  propagation  of  energy  along 
the  equator  in  the  late  summer  and  during  the  northeast  monsoon.  The  model 
produces  a  reversing  Somali  Current  with  southward  flow  during  the  northeast 
Monsoon  and  northward  flow  during  the  northeast  monsoon.  The  eddies  associated 
with  the  latter,  The  Southern  Gyre,  the  Great  Whirl  and  the  Socotra  eddy  are  all 
reproduced  in  the  solution.  The  northward  migration  of  the  southern  gyre  in  the 
final  phase  of  the  Monsoon  is  absent  with  the  wind  forcing  used  here. 

Considering  the  low  vertical  resolution  in  the  model,  the  observed  seasonal 
behaviour  of  the  undercurrents  in  the  Somali  Current  system  is  remarkably  well 
reproduced.  Linear  theories,  e.g.  Ligthill,  (1969),  McCreary  and  Kundu,  (19S5), 
show  that  the  first  few  vertical  modes  dominate  the  oceanic  response  to  the  wind, 
so  this  result  could  be  anticipated.  The  vertical  structure  observed  along  the  Somali 
coast  is  also  dominated  by  the  lowest  modes,  see  for  instance  Leetmaa  et  ai,  (1982), 
while  higher  modes  also  play  a  role  in  the  equatorial  flow  (Luyten  and  Swallow, 
1976). 

The  disappearence  of  the  Somali  coastal  undercurrent  in  the  spring  is 
found  to  be  caused  in  part  by  the  arrival  of  the  deep  equatorial  undercurrents  due 
to  semi-annual  zonal  wind  forcing.  The  regeneration  of  the  undercurrent  in  the  fall, 
which  agrees  with  the  observations,  is  associated  with  the  decay  of  the  Great  Whirl 
in  our  model  solution.  In  addition  to  these  new  results,  the  model  solution  shows  a 
complicated  pattern  of  eddies  between  the  equator  and  the  southern  Arabian  Coast 
in  the  deeper  layers  where  observations  are  sparse.  Since  there  is  no  thermodynamic 
forcing,  our  model  results  suggest  these  undercurrents  to  be  mainly  wind-driven. 
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The  generation  and  decay  of  the  Great  Whirl  is  described  in  detail  for 
year  21  of  the  model  solution.  During  the  initial  phase  of  the  summer  monsoon 
the  numerical  solution  show  the  characteristics  of  an  analytical  linear  vertical  mode 
solution  of  McCreary  and  Kundu,  (1985).  It  is  demonstrated  that  horizontal  shear 
instability  is  the  most  likely  mechanism  involved  in  its  formation.  This  corroborates 
earlier  findings,  e.g.,  Cox  (1979).  Maps  of  potential  vorticity,  relative  vorticity, 
kinetic  energy  and  available  potential  energy  are  presented  to  help  understanding 
the  dynamics  in  the  region. 

Many  physical  processes  are  involved  during  the  decay  of  the  Great  Whirl. 
The  generation  of  a  cyclonic  eddy  between  the  whirl  and  Socotra  is  followed  by 
southeastward  mutual  advection  of  the  eddies.  A  rapid  decrease  in  available  po¬ 
tential  energy  of  the  Great  Whirl  and  an  increase  of  cyclonic  kinetic  energy  in  the 
lower  layers,  suggest  that  baroclinic  instability  is  partly  responsible  for  the  decay 
in  the  3.5  layer  case. 

Two  experiments  where  the  duration  of  the  summer  monsoon  was  extended 
by  keeping  the  wind  stress  constant  after  July  16  and  August  16,  respectively, 
showed  a  similar  solution.  This  demonstrates  that  the  formation  of  the  cyclonic 
eddies  are  not  caused  by  wind  forcing,  but  is  due  to  redistribution  of  energy  in  the 
flow  itself. 

With  a  single  layer,  vertical  energy  exchange  is  impossible,  and  the  tran¬ 
sition  to  the  winter  monsoon  is  therefore  different.  Here  energy  is  transfered  from 
the  Great  Whirl  eastward  to  the  Socotra  eddy.  As  a  result  of  the  changed  decay 
mechanism  this  case  is  different  from  the  3.5  layer  solution  during  the  following 
winter  monsoon. 


5.  DISCUSSION  AND  CONCLUSIONS 

Many  authors  have  been  successful  in  modelling  some  aspects  of  the  re¬ 
sponse  of  the  Somali  Current  to  the  summer  monsoon,  but  no  model  has  been 
presented  to  date  that  is  capable  of  reproducing  the  observed  features  of  the  un¬ 
dercurrents.  It  is  clear  from  earlier  studies  that  realistic  coastlines  and  winds  are 
nessessary  to  obtain  the  observed  configuration  in  numerical  models,  e.g.,  Knox  and 
Anderson,  (1985).  Because  of  the  small  horizontal  scales  along  the  western  bound¬ 
ary,  the  internal  radius  of  deformation  must  also  be  resolved.  These  requirements 
made  it  nessessary  for  us  to  use  the  latest  supercomputer  technology  and  a  very 
efficient  code. 

The  fine  horizontal  resolution  enables  our  model  to  reproduce  more  of  the 
detailed  structure  in  the  Somali  Current  system  than  in  earlier  studies.  The  solu¬ 
tion  includes  the  coastal  upwelling  in  the  spring  and  generation  of  the  Great  Whirl 
at  the  correct  latitude  with  following  northward  migration.  This  has  previously 
been  modelled  for  the  upper  layer  using  reduced-gravity  models  with  fine  horizon¬ 
tal  resolution,  but  our  model  simulates  in  addition  the  seasonal  variability  of  the 
undercurrents  in  very  good  agreement  with  the  sparse  observations.  The  results 
are  surprisingly  good  in  spite  of  the  very  limited  vertical  resolution,  lack  of  bottom 
topography  and  thermodynamic  forcing. 

The  most  important  physical  process  in  this  model  compared  to  single¬ 
layer  simulations  is  the  vertical  transfer  of  energy,  which  allows  baroclinic  instability. 
During  the  transition  from  summer  to  winter  monsoon  the  results  from  the  model 
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suggest  that  several  mechanisms  play  a  role  in  the  decay  of  the  Great  Whirl.  These 
are  migration  through  eddy-eddy  interaction,  eddy-wave  interaction  and  downward 
energy  propagation,  indicating  that  baroclinic  instability  is  an  important  part. 

Compared  to  large  ocean  general-circulation  models,  the  advantage  of  our 
simple  model  is  that  it  requires  relatively  modest  computer  resources.  This  allows 
experiments  with  additional  physics,  for  instance  inclusion  of  a  mixed  layer  or  in¬ 
terfacial  friction  between  layers,  or  simulations  with  other  wind  fields  or  geometries. 
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Appendix  A:  Derivation  of  model  equations 

In  this  section  the  derivation  of  the  transport  equation  will  be  presented 

in  more  detail.  Assuming  that  the  aspect  ratio  for  the  motion  is  small,  i.e.,  that 

the  horizontal  scale  of  motion  is  much  larger  than  the  vertical,  we  can  neglect 

the  vertical  velocity  everywhere,  except  in  the  continuity  equation.  The  governing 

equations  in  spherical  coordinates,  which  can  be  found  in  textbooks  (for  instance 

Pedlosky,  1979  or  Gill,  1982)  can  in  that  case  be  written  as  follows.  The  velocity 

component  in  the  east-west  direction,  positive  towards  east,  must  satisfy 


Du  /  u  \  .  1  dp  d  ( 4  <9u\ 

+  ^ )t,smS -  +  ai  (Aval)  +  5 


(A  —  1) 


where  the  total  (or  substantial)  derivative  is  defined  as 


D  &  _  „ 

—  =  —  +  v  •  V 
Dt  dt 


(A  -  2) 


Here  v  is  the  horizontal  velocity,  Q  the  rate  of  rotation  of  the  Earth  and  Ay  a 
vertical  kinematic  viscosity  coefficient.  Assuming  a  constant  eddy  viscosity  Afj, 
the  horizontal  friction  term  for  the  u  equation  is  given  by 

G°  =  Ah  V2u - o  ^  Ju(l  -  2cos2#)  +  2sin#  — -  (A  —  3) 

crcos^y  l  ,  d<p 


This  form  ensures  that  a  fluid  in  solid-body  rotation  around  an  axis  through  a  sphere 
does  not  experience  any  horizontal  shear  stress  in  spite  of  the  fact  that  the  velocity 
varies  around  the  sphere.  For  instance  a  zonal  flow  with  a  velocity  distribution  with 
latitude  given  by  u  =  uocos^  should  have  no  friction.  The  terms  in  addition  to  the 
Laplacian,  which  often  is  seen  used  alone  in  many  models,  are  needed  to  satisfy  this 
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10 

condition.  It  is  easily  seen  from  a  scale  analysis,  that  these  terms  are  very  small 
except  near  the  poles.  However,  for  the  sake  of  generality  they  are  retained  in  this 
model. 


For  the  velocity  component  in  the  north-south  direction,!),  we  have 


dv  f  u  \  ■  n  1  dp  d  ( 

Tt  +  (2n  +  ^ H 6  =  +  s(Av 


<9u\ 
dz  ) 


+  S6 


with  the  horizontal  friction  term 


gb  =  ah 


v-777e[v{1-2cos2e)  +  2siDe^. 


In  each  layer,  we  will  assume  that  the  density  is  constant,  i.e., 


(A -4) 


(A -5) 


^  =  0 
Dt 


(A -6) 


This  implies  that  the  conservation  of  mass  is  governed  by  the  continuity  equation 

(A -7) 


1  du  Id,  ..  dw  n 

- o'ol  + - a  ^E(vcosd)  +  ~  0 

acosd  o<j)  acosd  dd  dz 


where  w  is  the  vertical  velocity  component. 

To  derive  the  transport  equations  we  integrate  across  a  layer  of  constant 
density  and  velocity.  For  instance  we  define  the  zonal  volume  transport 


over  a  layer  of  thickness 


ZT 

U  =  J  udz 

zb 


H  =  zj  -  zb 


(A -7) 


(A -8) 


The  meridional  transport  component  V  is  defined  similarly. 

Let  us  first  consider  the  continuity  equation.  The  first  term  becomes,  using 
Leibnitz  rule 


I 


1  du ^ 
acosd  d<f> 


zb 


1 

acosd 


dU 

d<p 


dH 
U  d4> . 


(A -9) 


107 


The  second  term  yields 


fid.  m  .  f  f  tan#  1  dv\ 

/  - 7— (ncosWIa^  =  / - u  H — —  <22 

J  acosO  d0'  J  \  a  a  89  J 

ZB  Zb 

tan 0  1  dV  v  dH 

a  +  a  dO  a  dO 

Integrating  the  third  term  eliminates  w  from  the  system  of  equations: 

Zt 


(A  -  10) 


I 


ZB 


^ dz  =  w(zT)  -  w{zB) 


Dzj'  Dzb  DH 

“*  ~dT 


Dt  Dt  Dt 

where  the  last  equality  is  possible  due  to  the  lack  of  vertical  shear  in  the  horizontal 
velocities.  Adding  the  right  hand  sides  of  equations  (A  —  9)  to  (A  —  11)  gives  the 
continuity  equation  (2  —  7),  viz. 


dHi 


+ 


1 


dUj  d  ...  ■ 

1 1*  +  m{  icosB) 


dt  acosO 
The  acceleration  term  in  the  u  equation  is 


0 


(A  -  12) 


Du  du 


Dt  dt  a  cos  0  d<j>  a  dO 


u  du  v  dv 

d-  r\ 


(A  -  13) 


The  first  term  becomes 

ZT 


/ 

ZB 


du  ,  dU  dH 
Ttdz  =  -dl~u~m 


(A -14) 


the  third  can  be  written,  using  v  constant  over  the  layer: 

Zt 


v  f  du  1  ' 

ZJ  MdZ  =  -al 


1  r  dvU  ..dv  dH  1 

W-vm~mTe 


ZB 


lr d ,uvs  dv 


(A -15) 
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From  this,  it  is  easily  seen  that  the  second  term  becomes 

ZT 


/ 


2B 


a  cos  9  d<f>'~  acos0ld4>  H  d<f> 


udu  1 

■dz  = 


Using  the  continuity  equation  (A  —  12),  we  find 


-T 

/ 


ZB 


Du  dU  1  d  U2  1  d  UV  UVUnd 
Dt Z  ~  dt  +  acosedctD  H  ]  +  ad9{  H  )  aH 


Defining  the  Coriolis  parameter 


(A -16) 


(A -17) 


f  =  20.  sin  9 


(A -18) 


we  find  that  the  total  planetary  acceleration  is 


7(  v  + 

J  \  a  cos  9  / 


UV  tan  9 
aH 


ZB 


(A -18) 


which  completes  the  left  hand  side  of  equation  (2  —  2). 

The  horizontal  friction  term  (2  —  3)  is  easily  obtained  from  (A  —  3)  by 
using  the  definitions  u  =  U/H  and  v  =  V/H ,  and  multiplying  by  H ,  again  using 
the  assumption  that  we  do  not  have  any  vertical  variation  of  the  horizontal  fields 
over  a  layer.  The  vertical  friction  is  integrated  as 


ZT 


1 


ZB 


ddu  (du  du 

7~d2— O0  =  djl  —  —  7T- 

OZ  OZ  \uz  ZT  OZ  ZB 


)  =  (tt  -  tb)/p 


(A -19) 


The  terms  in  the  prediction  equation  (2  —  4)  for  V  are  obtained  in  the  same  fashion. 

We  have  not  yet  considered  the  pressure  gradients.  From  figure  1,  which 
represents  the  case  of  N  =  4,  we  see  that  the  hydrostatic  pressure  in  the  upper 
layer  must  be  given  by 


Pl{z)  =pa~  gp\z 


(A -20) 
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recalling  that  pa  is  the  atmospheric  pressure  and  the  z  coordinate  is  negative  in  the 
ocean  and  zero  at  the  sea  surface  rj{(f>^9).  For  deeper  layers  the  top  is  at 


j- 1 

vj  =  -E^' 

i=l 

so  we  find 


Pj{z)  =  pj_i(zTj)  -  gpj(z  -  zTj) 


(A -21) 


(A -22) 


where  the  pressure  at  the  top  of  layer  j  equals  the  pressure  at  the  bottom  of  layer 
j  ~  1,  e.g., 

j- 1 

Pj-\(2Tj)  =  pa  +  ^gpiHi  (A -23) 

i=l 

We  need  to  find  the  horizontal  pressure  gradient,  vertical  integrated  over 
each  layer.  It  is  straight  forward,  but  requires  a  little  algebra,  to  integrate  first  and 
then  use  Leibnitz  formula  to  obtain  this  quantity.  However,  if  we  keep  in  mind  that 
the  vertical  coordinate  varies  in  the  horizontal  direction  so  that 

N 

Wz  =  -Vp  =  -V  ^{D  +  Hi)  (A  -  24) 

i=l 


we  find,  taking  the  gradient  of  (A  —  22)  and  inserting  the  definitions  (A  —  21)  and 
(A -23): 

j- 1 

Vpj(z)  =  V(p„  -  E s(Pj  -  PiW  +  V7  (A  -  25) 

!=1 

which  integrates  vertically  by  multiplication  with  Hj  and  gives  the  pressure  terms 
as  written  in  (2  —  2)  and  (2  —  4). 


Appendix  B:  Open  boundary  conditions 

Most  oceanographic  models  cover  a  limited  area  of  the  world  ocean.  It  is 

therefore  nessessary  to  implement  boundary  conditions  across  parts  of  the  ocean. 
Early  models  used  walls  with  a  free-slip  condition  along  a  latitude  where  the  wind- 
stress  curl  is  assumed  to  be  zero.  For  linear  steady  models  this  is  a  natural  boundary, 
since  no  meridional  transport  occurs.  In  time-dependent  simulations  a  sponge  layer 
where  incoming  waves  are  dispersed  due  to  large  amount  of  friction,  has  often  been 
used.  Here  we  adopt  the  definition  of  an  open  boundary  by  Rped  and  Cooper, 
(1986),  to  be  a  boundary  where  a  wave  can  leave  the  computational  domain,  with¬ 
out  reflection,  and  no  information  from  the  outside  is  prescribed.  In  that  sense  the 
boundary  condition  described  below  is  not  strictly  open,  since  a  small  mount  of 
diffusion  and  relaxation  towards  a  far  field  has  been  added. 

First  the  calculation  of  the  vertical  modes  is  considered.  The  approach  is 
that  of  Lighthill  (1969)  except  we  are  interested  in  the  modes  for  the  case  where  the 
N-th  layer  is  at  rest.  For  simplicity  we  use  the  linearized  equations  on  an  f-plane 
to  derive  the  normal  modes. 

The  inviscid  equations  are  in  transport  form 


dUi 

1  dPj 

II 

1 

■'  -w 

CO 

3  4.  F- 
Pj  dx  3 

(B  —  1) 

dVj 

1  dPi  „ 

pi  +  Gj 

Pj  dy 

(B  —  2) 

where  the  vertical  integrated  pressure  is  given  by 

N- 1  _  j-1 

Pi  =  9«0J  (pj  Y,  eE7TiHi  -  B'i  -  «>"■•)  (B  -  3) 

i=l  FN  i=  1 
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and  Fj  and  Gj  is  the  external  forcing  on  layer  j.  The  continuity  equation  is 

(B  —  4) 


dHj  dUj  dVj 


dt  dx  dy 

We  proceed  by  inserting  the  expression  for  the  pressure  P,  (B  —  3)  into 
(B  —  1)  and  (B  —  2)  and  take  d  /dt.  Using  continuity  (B  —  4)  to  eliminate  dHi/dt 
from  these  equations  we  obtain 

AT?-  N~l  ,P)2tt.  P)2\a  , 

(B  —  5) 


and 


d2Uj 

dVj 

dFj 

N-l 

r-->  j 

'd2Uj 

d2Vj 

dt1 

f  dt 

dt 

+  9  2 L,  aP\ 

i= 1 

v  dx2  ~*" 

J 

dxdy 

d2Vj 

dUj 

dGj 

N-l 

(  d2  Uj 

d2V j 

dt2 

+  f0t 

dt 

+  9  2_^  aji  1 

t=l 

\dxdy  + 

J 

dy 2 

(B  —  6) 


where  the  matrix  aji  is  given  by 

«*■«*«  (B-7) 

V  pj  pN  / 

which  is  of  the  same  form  as  the  equations  given  by  Lighthill,  except  for  the  coeffi¬ 
cients  in  the  matrix  (B  —  7). 

The  solution  to  ( B  —  5)  and  (B  —  6)  can  be  expressed  as  a  sum  of  normal 
modes,  i.e., 

JV-1 

(B  —  8) 


N-l 

Uj  =  £ 

k= 1 


where  is  the  j-th  component  of  the  fc-th  eigenvector  of  the  matrix  ( B  —  7),  and 
Uk  the  amplitude  of  the  k- th  mode  given  by 


Uk  =  bkjUj 


(B  —  9) 


is  the  product  of  the  vertical  transport  in  each  layer  with  the  fc-th  row  of  the  matrix 
bkj,  which  is  the  inverse  of  the  matrix  akj ,  which  consists  of  the  N  —  1  eigenvectors 


of  at.;  as  columns. 
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The  N  —  1  eigenvalues  h k  are  the  equivalent  depths,  which  give  the  internal 
gravity  phase  speeds  of  the  vertical  modes  as 

4  =  9hk  (B  -  10) 

As  a  orthogonality  condition  we  will  require  that 

N- 1 

£  =  S„m  (B  —  11) 

j-  1 

where  8nrn  is  the  Kronecker  delta.  The  only  external  forcing  considered  in  this 
study  is  the  wind  stress  acting  as  a  body  force  across  the  upper  layer,  he., 

(Fi,Gl)  =  (i’,T*)ln  (B  -  12) 


which  implies  that  the  wind  stress  is  projected  onto  the  n-th  vertical  mode  as 

Q1  (tX,  ry)l  H\ 

E1 

j=l 

However  for  the  open  boundary  conditions  as  implemented  in  this  study,  the  scaling 
is  unimportant,  as  will  be  seen  in  what  follows. 

The  finite-difference  form  of  the  equations  (2  —  2)  and  (2  —  4)  cannot  be 
computed  along  open  boundaries,  since  the  values  nessessary  for  such  a  computa¬ 
tion  are  not  available.  To  some  extent  one-sided  difference  schemes  can  be  used, 
which  was  tried  in  some  tests,  but  it  is  not  possible  for  all  terms  in  the  equations. 
A  popular  method  is  to  use  an  extrapolation  of  interior  values  towards  the  bound¬ 
ary.  Here  a  Sommerfeld  radiation  condition  is  applied.  This  method  assumes  that 
disturbances  travel  as  free  waves  towards  the  boundary,  and  separates  the  forcing 
and  other  physics.  For  a  discussion  of  this  see  Rped  and  Smedstad  (1984)  and  the 
review  papers  by  Rped  and  Cooper  (1986,  1987).  We  illustrate  an  open  boundary 
to  the  east  with  waves  propagating  towards  the  boundary  in  the  x-  direction,  e.g., 


(B  —  13) 


(Fn,  Gn) 


3$  3$ 

~dt  +  3x  =  ° 


(B  —  14) 
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where  <$  can  be  any  of  the  dependent  variables,  and  c$  is  a  phase  speed. 

The  equation  above  can  be  solved  locally  for  c ^  near  the  boundary  by 
using  finite  differences.  If  we  have  the  solution  at  time  level  n,  we  compute  the 
phase  speed 

c$=~—  f",1 - (B  —  15) 

A|  nn'i  -  *b-2 

where  m  =  n  —  2  and  At*  =  2Af  when  the  leap-frog  scheme  is  used  while  m  =  n  —  1 
and  At*  =  At  for  the  Euler  forward  scheme.  Here  the  subscript  B  —  1  denoted  the 
variable  one  grid  point  from  the  boundary. 

The  new  value  at  the  boundary,  can  the  found,  for  instance  using 

the  modified  Orlanski  condition  suggested  by  Miller  and  Thorpe,  (1981): 

^5  c*  <0 

(1  —  c<j>)^!b  +  1  C<J>  >  0 

where  is  the  non-dimensionalized  phase  speed, 


*n+l  _ 


c$  =  min(c$— ,1)  (B  —  17) 

The  Camerlengo  and  O’Brien,(1980)  scheme  corresponds  to  always  using  the  max¬ 
imal  phase  speed,  he.,  =  1,  when  the  wave  motion  is  towards  the  boundary.  It 
should  be  noted  ,  that  the  original  schemes  by  Orlanski  and  by  Camerlengo  and 
O’Brien  prescribed  =  ♦S"1  when  c<|>  <  0.  The  tests  described  below  showed 

no  dependence  in  the  solutions  on  which  time  level  was  used. 

Associated  with  the  computation  of  the  boundary  value  $£j+1  it  is  common 
to  apply  a  smoothing  operation  as  done  below.  This  will  help  absorb  reflecting 
disturbances  and  reduce  small  scale  noise. 

The  set  of  equations  (B  —  15)  to  (B  —  17)  can  be  applied  to  any  of  the 
variables  U,  V  and  H ,  or  their  normal  mode  amplitudes  U,  V  and  H . 

The  following  test  case  was  considered  to  evaluate  the  boundary  conditions 
with  the  least  amount  of  reflection:  A  rectangular  equatorial  2.5  layer  ocean  with 
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open  southern  and  eastern  boundaries  was  forced  with  a  zonal  wind  stress  in  the 
western  part  of  the  basin  to  create  an  equatorial  Kelvin  wave.  The  densities  of  the 
three  layers  are  1025  kg/m^,  1028.2  kg/m^  and  1029  kg/m^  while  the  initial  layer 
thicknesses  of  the  two  active  layers  are  200  m  and  300  m,  respectively.  The  model 
equations  and  other  parameters  are  the  same  as  used  in  the  model  simulations  of 
the  main  paper,  except  that  the  basin  only  was  extended  to  80°  E. 

Figure  51  shows  the  solution  in  the  two  layers  after  10,  20  and  40  days.  The 
first  baroclinic  Kelvin  mode  is  passing  through  the  eastern  boundary  from  day  14 
to  day  30.  Since  the  second  vertical  mode  moves  much  slower  we  have  a  separation 
of  the  two  Kelvin  modes.  At  day  40  the  second  baroclinic  mode  is  moving  out  of 
the  domain.  We  see  reflections  into  weak  westward  propagating  Rossby  waves,  but 
no  coastal  Kelvin  are  propagating  poleward  along  the  western  boundary. 

Several  possible  combinations  of  the  equations  above  were  considered.  One 
group  of  experiments  applied  (B  —  15)  to  (B  —  17)  to  the  dependent  variables  layer 
by  layer.  At  the  eastern  boundary,  U  was  always  calculated  using  an  open  boundary 
condition.  Cases  where  H  and  V  was  calculated  using  the  actual  equations,  with 
onenotrh-south-sided  finite  differences  when  nessessary,  were  considered  in  addition 
to  cases  where  the  open  boundary  condition  were  used  for  these  variables.  Both  the 
modified  Orlanski  condition  (B  —  16)  and  the  Camerlengo-O’Brien  condition  were 
tested.  A  second  group  of  experiments  were  made  where  the  variables  considered 
were  transformed  into  normal  vertical  modes  prior  to  the  application  of  the  open 
boundary  condition. 

The  condition  which  show  the  minimum  reflection  of  the  equatorial  Kelvin 
wave  into  equatorial  Rossby  waves  was  the  Camerlengo-O’Brien  condition  when 
it  was  applied  to  U  and  V  and  H  was  computed  from  the  continuity  equation. 
However,  the  differences  found  between  the  different  cases  were  not  very  large,  and 


Figure  51.  Equatorial  Kelvin  wave  in  2.5  layer  model.  The  initial  depth  to  the 
bottom  of  layer  1  and  2  is  200  m  and  500m.  Top  to  bottom:  Day  10,  20  and 
40,  in  layer  1  (left)  and  layer  2  (right).  Eastern  and  southern  boundaries  are 
open. 
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all  cases  produced  reflections,  which  made  it  nessessary  to  apply  some  smoothing. 
A  two-dimensional  Hanning  filter,  i.e., 

®i,j  =  4  ®ij  +  1  + 

+  +  ^i+lj-l  +  $t-l,j+l  +  ^i+lj-t-l)  {B  -  18) 

was  applied  to  the  two  rows  or  columns  nearest  to  the  boundary  for  U,  V  and  H. 
For  the  boundary  points  a  one-sided  smoothing  was  applied  perpendicular  to  the 
boundary.  If  done  every  time  step  the  absorbtion  is  very  large,  in  fact  the  effective 
horizontal  viscosity  along  the  boundaries  is  two  orders  of  magnitude  higher  that 
in  the  interior,  and  the  choice  of  open  boundary  scheme  is  not  important.  With 
smoothing  less  frequent  than  every  10  time  step  reflections  starts  to  become  visible. 
As  a  compromise  the  Hanning  filter  was  applied  each  7^  time  step  in  this  study. 

In  the  Indian  Ocean  simulations  it  became  clear  that  the  open  boundary 
conditions  needed  further  modification  where  the  future  value  of  a  variable  depends 
on  the  conditions  outside  the  domain.  Along  the  east  coast  of  Africa,  coastal  Kelvin 
waves  travel  towards  the  equator,  and  the  same  type  of  waves  propagate  counter 
clockwise  around  Madagascar.  This  implies  that  the  layer  thicknesses  Hj  is  deter¬ 
mined  from  the  conditions  outside  the  domain.  With  this  modification  the  model 
failed  to  produce  outflow  along  these  coasts  as  seen  in  observations. 

A  relaxation  towards  observations  was  used  by  Cho  and  Clark  (1981). 
When  the  flow  was  directed  into  the  interior  of  their  domain,  they  adjusted  the 
solution  obtained  from  the  Orlanski  scheme  towards  the  far  field  conditions  outside. 
For  the  two  southermost  Hj  points  along  the  coast  of  East  Africa  and  Madagascar, 
we  apply  the  relaxation,  regardless  of  the  direction  of  the  flow, 


Hj  =  (1  —  e)Hj  -f-  cHftj 


(B  —  19) 
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where  H  ■  is  the  adjusted  value  and  H$j  is  the  initial  layer  thickness.  According 
to  the  atlas  by  Wyrtki  (1971),  Hqj  is  a  reasonable  choice  for  average  thicknesses 
in  the  southern  part  of  the  Indian  Ocean.  It  was  attempted  to  adjust  the  layer 
thicknesses  on  the  eastern  side  of  Madagascar  to  those  on  the  western  side  since 
pressure  differences  can  propagate  as  coastal  Kelvin  waves.  It  worked  well  for  the 
1.5  layer  model,  but  the  3.5  layer  case  developed  a  strong  seasonal  variation  in  the 
outflow  along  the  east  coast,  and  negatively  correlated  with  that,  in  the  flow  north 
of  Madagascar,  which  again  influenced  the  conditions  on  the  western  side  of  the 
island.  Since  the  amplitude  of  this  seasonal  variation  in  transport  increased  in  time 
and  may  be  due  to  some  artificial  resonance,  this  approach  was  not  used. 

The  relaxation  towards  the  far  field,  even  if  it  is  done  in  a  few  points  as 
in  this  study  does  put  restrictions  on  the  flow  condition.  However  some  constraints 
may  be  nessessary.  Hurlburt  and  Thompson  (1980)  found  that  recirculation  through 
an  open  boundary,  if  unrestricted,  can  attain  large  unrealistic  transports.  To  avoid 
this  problem,  Thompson  and  Schmitz  (1989)  who  used  a  modified  Orlanski  open 
boundary  condition,  imposed  a  damping  similar  to  (B  —  19)  at  all  inflow  points, 
and  scaled  the  total  outflow  uniformly  to  match  the  total  inflow. 

The  formulation  of  the  open  boundaries  does  not  constrain  the  inflow  and 
outflow  to  conserve  mass.  The  mass  loss  through  the  open  boundaries  correspond  to 
a  few  meters  of  layer  thickness  per  year  of  integration.  While  this  is  acceptable  for 
shorter  integrations,  a  compensation  was  needed  in  this  study.  In  each  time  step  the 
global  average  layer  thickness  anomaly  is  calculated  for  each  layer,  and  subtracted 
from  the  thickness  at  each  grid  point,  (Luther,  personal  communication).  This  can 
be  thought  of  as  a  vertical  diffussion  of  mass  through  the  layers.  Since  the  pressure 
gradients  are  unchanged  by  this  adjustment,  the  geostrophic  balance  is  not  affected. 
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It  should  be  clear  from  the  discussion  above,  that  the  problem  of  imple¬ 
menting  open  boundary  conditions  for  realistic  ocean  models  has  not  been  solved 
satisfactorily  at  this  time,  and  there  is  a  need  for  further  work  on  this  topic. 


